{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loopy: Dealing with Intermediate Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyopencl as cl\n",
    "import pyopencl.array\n",
    "import pyopencl.clmath\n",
    "import pyopencl.clrandom\n",
    "import loopy as lp\n",
    "\n",
    "from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Choose platform:\n",
      "[0] <pyopencl.Platform 'Portable Computing Language' at 0x7f79b0a7e6e8>\n",
      "[1] <pyopencl.Platform 'Intel(R) OpenCL' at 0x1af0e48>\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Choice [0]: \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set the environment variable PYOPENCL_CTX='' to avoid being asked again.\n"
     ]
    }
   ],
   "source": [
    "ctx = cl.create_some_context(interactive=True)\n",
    "queue = cl.CommandQueue(ctx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "grid = np.linspace(0, 2*np.pi, 2048, endpoint=False)\n",
    "h = grid[1] - grid[0]\n",
    "u = cl.clmath.sin(cl.array.to_device(queue, grid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two kernels: Finite Differences and a Flux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will implement computing the ODE right-hand side for Burgers' equation:\n",
    "$$\n",
    "\\frac{\\partial u}{\\partial t} + \\frac{\\partial }{\\partial x} \\left( \\frac{u^2}2 \\right) = 0,\n",
    "$$\n",
    "\n",
    "Now, it is likely that the code fore the derivative is initially separate from the code for the flux function $f(u):=u^2/2$.\n",
    "\n",
    "For simplicity, we will use central finite differences to build a kernel `fin_diff_knl` to take the derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "f: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [i] : 0 < i <= n }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "i: None\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for i\n",
      "  \u001b[36mout[i]\u001b[0m = \u001b[35m((-1)*(f[i + 1] + (-1)*f[i + -1])) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end i\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "fin_diff_knl = lp.make_kernel(\n",
    "    \"{[i]: 1<=i<=n}\",\n",
    "    \"out[i] = -(f[i+1] - f[i-1])/(2*h)\",\n",
    "    [lp.GlobalArg(\"out\", shape=\"2+n\"), ...])\n",
    "print(fin_diff_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------\n",
    "Next, make the flux kernel as `flux_knl`.\n",
    "* Use `j` as the loop variable.\n",
    "* Make sure to declare `f` and `u` to be the right size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "f: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [j] : 0 < j <= n }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "j: None\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for j\n",
      "  \u001b[36mf[j]\u001b[0m = \u001b[35mu[j]**2 / 2\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end j\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "flux_knl = lp.make_kernel(\n",
    "    \"{[j]: 1<=j<=n}\",\n",
    "    \"f[j] = u[j]**2/2\",\n",
    "    [\n",
    "    lp.GlobalArg(\"f\", shape=\"2+n\"),\n",
    "    lp.GlobalArg(\"u\", shape=\"2+n\"),\n",
    "    ])\n",
    "print(flux_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fuse the Kernels\n",
    "\n",
    "Next, fuse the two kernels together as `fused_knl`, using `lp.fuse_kernels([knl_a, knl_b])`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel_and_loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "f: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [i] : 0 < i <= n }\n",
      "[n] -> { [j] : 0 < j <= n }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "i: None\n",
      "j: None\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for i\n",
      "  \u001b[36mout[i]\u001b[0m = \u001b[35m((-1)*(f[i + 1] + (-1)*f[i + -1])) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end i\n",
      "for j\n",
      "  \u001b[36mf[j]\u001b[0m = \u001b[35mu[j]**2 / 2\u001b[0m  {id=\u001b[32minsn_0\u001b[0m}\n",
      "end j\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "#clear\n",
    "fused_knl = lp.fuse_kernels([fin_diff_knl, flux_knl])\n",
    "\n",
    "print(fused_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the generated code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/andreas/src/loopy/loopy/kernel/creation.py:1803: LoopyWarning: in kernel loopy_kernel_and_loopy_kernel: The single-writer dependency heuristic added dependencies on instruction ID(s) 'insn_0' to instruction ID 'insn' after kernel creation is complete. This is deprecated and may stop working in the future. To fix this, ensure that instruction dependencies are added/resolved as soon as possible, ideally at kernel creation time. (add 'single_writer_after_creation' to silenced_warnings kernel attribute to disable)\n",
      "  warn_with_kernel(kernel, \"single_writer_after_creation\",\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine lid(N) ((int) get_local_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine gid(N) ((int) get_group_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mif __OPENCL_C_VERSION__ < 120\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mpragma OPENCL EXTENSION cl_khr_fp64: enable\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mendif\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\n",
      "__kernel \u001b[36mvoid\u001b[39;49;00m \u001b[32m__attribute__\u001b[39;49;00m ((reqd_work_group_size(\u001b[34m1\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m))) loopy_kernel_and_loopy_kernel(__global \u001b[36mdouble\u001b[39;49;00m *__restrict__ out, __global \u001b[36mdouble\u001b[39;49;00m *__restrict__ f, \u001b[36mfloat\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m h, \u001b[36mint\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m n, __global \u001b[36mdouble\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m *__restrict__ u)\n",
      "{\n",
      "  \u001b[34mfor\u001b[39;49;00m (\u001b[36mint\u001b[39;49;00m j = \u001b[34m1\u001b[39;49;00m; j <= n; ++j)\n",
      "    f[j] = (u[j] * u[j]) / \u001b[34m2.0\u001b[39;49;00m;\n",
      "  \u001b[34mfor\u001b[39;49;00m (\u001b[36mint\u001b[39;49;00m i = \u001b[34m1\u001b[39;49;00m; i <= n; ++i)\n",
      "    out[i] = (-\u001b[34m1.0\u001b[39;49;00m * (f[\u001b[34m1\u001b[39;49;00m + i] + -\u001b[34m1.0\u001b[39;49;00m * f[-\u001b[34m1\u001b[39;49;00m + i])) / (\u001b[34m2.0\u001b[39;49;00m * h);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fused_knl = lp.set_options(fused_knl, write_cl=True)\n",
    "evt, _ = fused_knl(queue, u=u, h=np.float32(1e-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, eliminate the separate loop to compute `f`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel_and_loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [i] : 0 < i <= n }\n",
      "[n] -> { [] : n > 0 }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "i: None\n",
      "---------------------------------------------------------------------------\n",
      "SUBSTITUTION RULES:\n",
      "f_subst(j) := u[j]**2 / 2\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for i\n",
      "  \u001b[36mout[i]\u001b[0m = \u001b[35m((-1)*(f_subst(i + 1) + (-1)*f_subst(i + -1))) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end i\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "fused_knl = lp.assignment_to_subst(fused_knl, \"f\")\n",
    "print(fused_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let's take a look at the generated code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine lid(N) ((int) get_local_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine gid(N) ((int) get_group_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mif __OPENCL_C_VERSION__ < 120\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mpragma OPENCL EXTENSION cl_khr_fp64: enable\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mendif\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\n",
      "__kernel \u001b[36mvoid\u001b[39;49;00m \u001b[32m__attribute__\u001b[39;49;00m ((reqd_work_group_size(\u001b[34m1\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m))) loopy_kernel_and_loopy_kernel(__global \u001b[36mdouble\u001b[39;49;00m *__restrict__ out, \u001b[36mfloat\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m h, \u001b[36mint\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m n, __global \u001b[36mdouble\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m *__restrict__ u)\n",
      "{\n",
      "  \u001b[34mfor\u001b[39;49;00m (\u001b[36mint\u001b[39;49;00m i = \u001b[34m1\u001b[39;49;00m; i <= n; ++i)\n",
      "    out[i] = (-\u001b[34m1.0\u001b[39;49;00m * ((u[\u001b[34m1\u001b[39;49;00m + i] * u[\u001b[34m1\u001b[39;49;00m + i]) / \u001b[34m2.0\u001b[39;49;00m + -\u001b[34m1.0\u001b[39;49;00m * ((u[-\u001b[34m1\u001b[39;49;00m + i] * u[-\u001b[34m1\u001b[39;49;00m + i]) / \u001b[34m2.0\u001b[39;49;00m))) / (\u001b[34m2.0\u001b[39;49;00m * h);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fused_knl = lp.set_options(fused_knl, write_cl=True)\n",
    "evt, _ = fused_knl(queue, u=u, h=np.float32(1e-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transform the kernel for execution\n",
    "\n",
    "For easier transformation, renumber the loop to start at 0, using `affine_map_inames(kernel, old_inames, new_inames, equations)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel_and_loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [inew] : 0 <= inew < n }\n",
      "[n] -> { [] : n > 0 }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "inew: None\n",
      "---------------------------------------------------------------------------\n",
      "SUBSTITUTION RULES:\n",
      "f_subst(j) := u[j]**2 / 2\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for inew\n",
      "  \u001b[36mout[1 + inew]\u001b[0m = \u001b[35m((-1)*(f_subst(1 + 1 + inew) + (-1)*f_subst(-1 + 1 + inew))) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end inew\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "fused0_knl = lp.affine_map_inames(fused_knl, \"i\", \"inew\", \"inew+1=i\")\n",
    "\n",
    "print(fused0_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, map the kernel to OpenCL axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel_and_loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [inew_outer, inew_inner] : inew_inner >= 0 and -128inew_outer <= inew_inner <= 127 and inew_inner < n - 128inew_outer }\n",
      "[n] -> { [] : n > 0 }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "inew_inner: l.0\n",
      "inew_outer: g.0\n",
      "---------------------------------------------------------------------------\n",
      "SUBSTITUTION RULES:\n",
      "f_subst(j) := u[j]**2 / 2\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "for inew_outer, inew_inner\n",
      "    \u001b[36mout[1 + inew_inner + inew_outer*128]\u001b[0m = \u001b[35m((-1)*(f_subst(1 + 1 + inew_inner + inew_outer*128) + (-1)*f_subst(-1 + 1 + inew_inner + inew_outer*128))) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "end inew_outer, inew_inner\n",
      "---------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "gpu_knl = lp.split_iname(fused0_knl, \"inew\", 128, outer_tag=\"g.0\", inner_tag=\"l.0\")\n",
    "print(gpu_knl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, precompute the fluxes locally in each workgroup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------------------------------------------------\n",
      "KERNEL: loopy_kernel_and_loopy_kernel\n",
      "---------------------------------------------------------------------------\n",
      "ARGUMENTS:\n",
      "h: ValueArg, type: <auto/runtime>\n",
      "n: ValueArg, type: <auto/runtime>\n",
      "out: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "u: type: <auto/runtime>, shape: (2 + n), dim_tags: (N0:stride:1) aspace: global\n",
      "---------------------------------------------------------------------------\n",
      "DOMAINS:\n",
      "[n] -> { [inew_outer, inew_inner, j_outer, j_inner] : inew_outer >= 0 and 0 <= inew_inner <= 127 and inew_inner < n - 128inew_outer and j_inner >= 0 and -128j_outer <= j_inner <= 129 - 128j_outer and j_inner <= 127 and j_inner <= 1 + n - 128inew_outer - 128j_outer }\n",
      "[n] -> { [] : n > 0 }\n",
      "---------------------------------------------------------------------------\n",
      "INAME IMPLEMENTATION TAGS:\n",
      "inew_inner: l.0\n",
      "inew_outer: g.0\n",
      "j_inner: l.0\n",
      "j_outer: None\n",
      "---------------------------------------------------------------------------\n",
      "TEMPORARIES:\n",
      "f_subst_0: type: <auto/runtime>, shape: (j:130), dim_tags: (N0:stride:1) scope:auto\n",
      "---------------------------------------------------------------------------\n",
      "INSTRUCTIONS:\n",
      "  for inew_outer, j_outer, j_inner\n",
      "↱       \u001b[36mf_subst_0[j_inner + j_outer*128]\u001b[0m = \u001b[35mu[128*inew_outer + j_inner + j_outer*128]**2 / 2\u001b[0m  {id=\u001b[32mf_subst\u001b[0m}\n",
      "│   end j_outer, j_inner\n",
      "│   for inew_inner\n",
      "└     \u001b[36mout[1 + inew_inner + inew_outer*128]\u001b[0m = \u001b[35m((-1)*(f_subst_0[2 + inew_inner] + (-1)*f_subst_0[inew_inner])) / (2*h)\u001b[0m  {id=\u001b[32minsn\u001b[0m}\n",
      "  end inew_outer, inew_inner\n",
      "---------------------------------------------------------------------------\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine lid(N) ((int) get_local_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mdefine gid(N) ((int) get_group_id(N))\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mif __OPENCL_C_VERSION__ < 120\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mpragma OPENCL EXTENSION cl_khr_fp64: enable\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\u001b[36m#\u001b[39;49;00m\u001b[36mendif\u001b[39;49;00m\u001b[36m\u001b[39;49;00m\n",
      "\n",
      "__kernel \u001b[36mvoid\u001b[39;49;00m \u001b[32m__attribute__\u001b[39;49;00m ((reqd_work_group_size(\u001b[34m128\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m, \u001b[34m1\u001b[39;49;00m))) loopy_kernel_and_loopy_kernel(__global \u001b[36mdouble\u001b[39;49;00m *__restrict__ out, \u001b[36mdouble\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m h, \u001b[36mint\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m n, __global \u001b[36mdouble\u001b[39;49;00m \u001b[34mconst\u001b[39;49;00m *__restrict__ u)\n",
      "{\n",
      "  __local \u001b[36mdouble\u001b[39;49;00m f_subst_0[\u001b[34m130\u001b[39;49;00m];\n",
      "\n",
      "  \u001b[34mif\u001b[39;49;00m (\u001b[34m1\u001b[39;49;00m + -\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + -\u001b[34m1\u001b[39;49;00m * lid(\u001b[34m0\u001b[39;49;00m) + n >= \u001b[34m0\u001b[39;49;00m)\n",
      "  {\n",
      "    f_subst_0[lid(\u001b[34m0\u001b[39;49;00m)] = (u[\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + lid(\u001b[34m0\u001b[39;49;00m)] * u[\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + lid(\u001b[34m0\u001b[39;49;00m)]) / \u001b[34m2.0\u001b[39;49;00m;\n",
      "    \u001b[34mif\u001b[39;49;00m (\u001b[34m1\u001b[39;49;00m + -\u001b[34m1\u001b[39;49;00m * lid(\u001b[34m0\u001b[39;49;00m) >= \u001b[34m0\u001b[39;49;00m && -\u001b[34m127\u001b[39;49;00m + -\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + -\u001b[34m1\u001b[39;49;00m * lid(\u001b[34m0\u001b[39;49;00m) + n >= \u001b[34m0\u001b[39;49;00m)\n",
      "      f_subst_0[\u001b[34m128\u001b[39;49;00m + lid(\u001b[34m0\u001b[39;49;00m)] = (u[\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + \u001b[34m128\u001b[39;49;00m + lid(\u001b[34m0\u001b[39;49;00m)] * u[\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + \u001b[34m128\u001b[39;49;00m + lid(\u001b[34m0\u001b[39;49;00m)]) / \u001b[34m2.0\u001b[39;49;00m;\n",
      "  }\n",
      "  barrier(CLK_LOCAL_MEM_FENCE) \u001b[37m/* for f_subst_0 (insn depends on f_subst) */\u001b[39;49;00m;\n",
      "  \u001b[34mif\u001b[39;49;00m (-\u001b[34m1\u001b[39;49;00m + -\u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + -\u001b[34m1\u001b[39;49;00m * lid(\u001b[34m0\u001b[39;49;00m) + n >= \u001b[34m0\u001b[39;49;00m)\n",
      "    out[\u001b[34m1\u001b[39;49;00m + \u001b[34m128\u001b[39;49;00m * gid(\u001b[34m0\u001b[39;49;00m) + lid(\u001b[34m0\u001b[39;49;00m)] = (-\u001b[34m1.0\u001b[39;49;00m * (f_subst_0[\u001b[34m2\u001b[39;49;00m + lid(\u001b[34m0\u001b[39;49;00m)] + -\u001b[34m1.0\u001b[39;49;00m * f_subst_0[lid(\u001b[34m0\u001b[39;49;00m)])) / (\u001b[34m2.0\u001b[39;49;00m * h);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "precomp_knl = lp.precompute(gpu_knl, \"f_subst\", \"inew_inner\", fetch_bounding_box=True, default_tag=\"l.auto\")\n",
    "print(precomp_knl)\n",
    "precomp_knl = lp.tag_inames(precomp_knl, {\"j_outer\": \"unr\"})\n",
    "precomp_knl = lp.set_options(precomp_knl, return_dict=True)\n",
    "evt, _ = precomp_knl(queue, u=u, h=h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the PDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3iUxdqH79lN7z0hIYWSBEJLIPTeUQRERIpdBBuCiIKCFJGiFGkKokewKyAiSAfpTQglBAihJqT3nmx2szvfH7smWeI5x/MZCNH3vq69kp3nLTNZeH8785QRUkoUFBQUFP65qGq7AwoKCgoKtYsiBAoKCgr/cBQhUFBQUPiHowiBgoKCwj8cRQgUFBQU/uFY1HYH/j94eHjIoKCg2u6GgoKCQp3izJkzWVJKzzvb66QQBAUFERUVVdvdUFBQUKhTCCES/qhdWRpSUFBQ+IejCIGCgoLCPxxFCBQUFBT+4ShCoKCgoPAPRxECBQUFhX84NSIEQoi1QogMIcTFf2MXQogVQojrQogLQojWVWxPCyGumV5P10R/FBQUFBT+PDU1I/gCGPAf7A8AwabXOGA1gBDCDZgFtAfaAbOEEK411CcFBQUFhT9BjeQRSCkPCyGC/sMhQ4CvpLHm9UkhhIsQoh7QA9grpcwBEELsxSgo39dEvxT+N8rK9dzKKiY+q5i8Eh0FGh0GCTYWKpztLPF3tSPQQY+nNhlyb0FxFpRrQK8FK0ewcQbXIAqcfLihzSGjJIOs0izK9GUAWKms8LD1wNvOi4bSE8vkDHRJSegLCjEUF4FKjcrGBrWbG5YB/mhd/CgqVVGcq6GkQIfBYEAawMpWja2jFY5u1qgtCynOzaIoO4uSgnykwYDBoMfazgF7FxdcvOvhERCI2sKylv+6Cgr3L/cqocwPSKzyPsnU9u/aqyGEGIdxNkFAQMDd6eU/jHK9gd9u5XD4aibHbmQRm1qI3mC+P4U1WnqooumuisZRdRVPVVK16+SoVByzs+GorS0XrK1Isqz+0PXOkUTclDRNlDgnSlKKze0SKHKoT7ZbGHkuwRQ66tFZZpkfI/XI8lT05fEYylOQ5RmA9r+OU6W2oF5wCI0jOxDSsStOHtUSKxUU/tHUmcxiKeWnwKcAkZGRym46f4H4rGK+OZnAlugUMgvLsFQLIgJcebF7Q0K8HWnkYY9P3hmcL32DxfXdCF0xeitH8twjOGM7mDOl3hzMsiZWJFLucgnsboGQuFk508YhkGFYEVKQhff1i9hd1qFNtEObozbe3McTTcdG3PC14pJ9PtGaYtwKwwkqaIOt1gkAVw9LAh21OJWlYnnrAvnXz5BuZyDT2Y5yASqVGvf6QTi6d0Zt5U1Rri056QIh7PBt7EZ4nwC8G9pQkp9HTnIiaTeukXDhHIe+Wcvh774gpH1nOjwyAo+AoNr7EBQU7iPulRAkA/5V3tc3tSVjXB6q2n7wHvXpH8fF5HxWHbzOzotpWKgEPUK9eCTCj+6hnthZWYC+HGI2wJYVkBkLNi7QagQ0HYw6qAvuaktkaRZHY78l/srXlOsKsVd5UpbXl/zsEFztG9OhezADDKkU/PoDhXtvUaTXY+PvhFfDDBx98rFq1YTyrq9gebseRYeScUguQqgh3yeZkzbbSHO7zpCWA+lefxhJJ7OJzsohz9cGK0sr/LR63OOT8ba2w/exZ3EeMhihMrq5ivPLiPstjYsHk9m55jK+wS50GR5MaMcgQjt2BSAvLZXofTuJ+XU3V08eo2WfAXQd/QzWdna1+bEoKNQ6oqa2qjT5CLZJKZv/gW0gMB54EKNjeIWUsp3JWXwG+D2K6CzQ5nefwb8jMjJSKrWG/jwpeaUs3HWFn8+n4GhjwZMdAnmmcxBejjbGAwwGuPQTHFwA2dfBuwW0fwGaDwMr40MyqzSLT6I/YfO1zegMOvoE9mFk6EgifSLRG2BXTCqHvv2FTkd/pmluAnoHRzwfexTXkSOxCgiA0lxKj37Fhf3xXMzrhkY64+FrTVi3QILbemNjb8ntgtus3bGEvOMX8cuwQSUF9cNa0LJXPxq374SFpRUlJ0+SuWw5pdHR2LZuje+C+VgFBlaM1aA3cPloCqe23aKsuJzIgUG0GRCISl0ZF1FaWMDxjd8RvWcHTp6eDJwwhXrBoff0M1FQqA2EEGeklJHV2mtCCIQQ32P8Zu8BpGOMBLIEkFJ+IoQQwEcYHcElwLNSyijTuc8B00yXmielXPff7qcIwZ+jXG/gsyO3WP7rVQwSxnZtwIvdG+FoU2UNP+0ibJ8MiSfBKwx6ToMmD4EQAJSWl/LVpa9Ye3EtWr2WIY2H8GzzZwl0qnz4lkRFkbF4CaXnz1Pu4cWmJn343qU5nZvXZ86QZnjaWHF+322i9yWiK9MTVC+bcN1H+NrdQjz4PrLlSG5Fn+HkT+tJvXoFaydHEgN1HHG7QbPGkSzougAPW4+K+0mDgfzNP5O+cCHodPjMnoXz4MFmY9cU6zj8w1WunU7HL9SVAeOaY2Nv7rtIjotlx8pFFOfm8OCrbxDSoctd+BQUFO4f7qoQ3GsUIfjvXM8oZPLGC0Qn5tG/mTczHgqjvmuVJRBtCRyYBydXg60L9J0DrUaDqvKb8+Gkw7x38j3SitPoHdCb11q/RpBzUIVdl5JCxuLFFOzYiYW3Nx4vvYjzI49gUFvwxfF4lu6+SiuNii46K2SZgUatPWk3qCFu9ewh5yb8/AqJsZc4UhhBaq4eRw9P2g0ZTvMefVBZWvDj1R9ZdHoRDlYOLOq2iEgf83+/utRUkt98k9KoM7iNeQ6vyZMrlop+J/Z4Kge/vYKThy2DXm2Fk4etmb20qJCfP5hDyrUr9Bv3Ki169au5D0FB4T5DEYJ/EBtOJ/LOlovYWamZM6Q5g1rWQ5i+4QOQGg2bxkJWHLR5BnrPAju3CnNWaRYLTy1kZ/xOGjk3YkbHGbTxblNhlzod2Z9/TtYna0BK3J9/Hvfnx6CyrXzIJlzK5sB3cRRna7hlocfQ3JlZT0fgbGv8Vp6dlMihrz/n1vkoHCzK6BhQTLPxq1H7tTIby9Xcq0w+OJnEwkRmdZzF0OChZnZZXk76/Pnkfvc9jv3747doIcLKyuyYlGt57Fh9AStbCx5+PQInd3Mx0GnL2LpkPgnR53ho0lRC2nf+f/3dFRTudxQh+Aeg0emZteUS66MS6dzYnaUjwiv9AGD0BRxfAfvngr0HPLwaGvU0u8b+2/uZdXwWxbpixrYcy/PNn8dSXbmkoomNJWXadMpiY3Hs1w/vqVOw9KuM+C3ILuXYxuvcPJ+Js5ctXR4LZld2Pot3x+HtZMOq4U3JOfIL53Zvw8rGlvZDHyO8uQ+Wm8dAaR4M+wyaDjLrU5G2iMmHJnM85Tjjw8czruU4M2GTUpKz7gsyFi7EoU9v6i9dirgjhDUjoYAty85jY2/BsCmR2DmZi4VOo2HjvHfIuHmdEe9+QL3Gis9A4e+HIgR/czILy3j+y9NEJ+UzvmdjJvUNQa2qMgvQFMDPL8GVbRA2BB5aZjYL0JRrWBy1mPVx62nq1pT3u75PQ5eGFXap1ZL1ySdkffoZahcXfGbNxKlv30q7QRJzKIkTm28AEPlgEOG9A1BbGpdqziTksHDlN7RIOYKtXkPLPv3pPOJJ7JycjRcoTIcfRkPKWXhoqXGmUgWdXses47P45eYvPNv8WSa1nmQ+ywFyvvmW9LlzjTODD5cg1Goze/qtAn7+8Cwe/o48PCmiom+/U1pUyDdvTURKyZPvL8fW0el/+xAUFO5z/p0QKEXn/gbcyCzikdXHuJpexKdPtuGN/qHmIpB1Df7VG+J2woAPYPiXZiJwK/8Wo7aPYn3cep4Ke4pvHvzGTAS0iYnEP/4EWatW4zzwQRpt+8VMBPIzS/l56TmOrL+Gb7Aro2d3oM2AoIoHbUFWBvHfLqNd4h609u6s9xtGQpMHK0UAwNEbnt4KjXrDLxPh8CKzMVqqLZnXZR4jQkew7uI6VkevrvZ3cHvicbymTqVw924yFi2uZvdu4ETvZ8JIu5nPwe/jqtltHRwZNOltSvJy2fvpR//9D6+g8DehziSUKfwxZxJyGfPladRC8MO4DrTydzE/4OZBWP8kqC3hqS3QoKuZ+XDSYaYenoqlypLVfVbTxc88ciZ/+3bSZs4CtRq/Fctx6lfpTJVSculwMsd+uoFKQM8nm9C0U6U/QkrJhX27OPztWqRB0uu5Fwnt3p+sTReYuz2W3BItb/QLrfxmb2UPo76HLa8Yl68Q0O2NivsJIZjWfhpl+jJWR6/G3tKep5uZ1yl0f/YZdMnJ5HzxBVZBQbiOHGFmb9zGi+zkIKJ2xOPfxJWQdj5mdp9GwXQc/jhHv/+Sa78dJ7h9pz/7USgo1FkUIajDnLqVwzPrTuHlaM2Xz7Uj0N3e/IALG43LQR7BMHoDuFTm9Ekp+SzmMz469xFN3JqwrOcyfB18K+wGjYa0uXPJ/3ETtuHh+C1ZbOYLKC3Ssv/LWOJjsvEPc6PnE01wdKv0R+Slp7H30xXcvniBgBbh9Bv3Ks5e3gCsHNUaZ9sYPj5wgyJNObMGNUP1+wxGbWn0XUgJ+98DCxvoNL7iuiqhYnbH2ZToSlgStYT6jvXpHdDbbNjeb01FezuBtHnzsAlrim3Llmb2tgODSLqSy6Hv4vBp6FwtkijyoaHEHT/Mr+s+IaBFuJJwpvD3R0pZ515t2rSR/3RO3MiSTWfslL0WH5Dp+aXVDzi2QspZTlKufUDKklwzk6ZcIycfnCybf9FcTjk0RZboSszs2pQUefORYfJyaBOZvuRDadBqzezJV3PkuqlH5apX9svo/belwWCosBn0enlmx1a57MlH5IqnH5XR+3aa2SuOMxjkvO2XZeDUbfLNjeelXn/HMeU6Kdc/aRxD1Lpq55fqSuXobaNl22/ayktZl6rZy/Py5NWePeW13n1keUFBNXt+ZolcM/Gg3Lr83B/2L+XqFbn4sYHy6PpvqtkUFOoqQJT8g2eq4iOog5y4kc2z607j52LL9+M64OVUJTJIStj3Lux5B8Iehid+MuYJmMgvy2fcnnHsjt/N621e5/2u72NrUfmNuCQqiluPDkcbH0/9Vavwen1SRQSOwSA5vf0WP394DgsrFY9OiaRlT/+KpZ3ivFw2LZjFgS/W4N+0OU8vXkXL3gOqOXXBuMzz9gNNmNCrMRuikpiz7TKyauCC2gKGfQ6N+8C21+H6PrPzbSxsWN5rOS7WLkzYP4FcTa6ZXe3sjN+SJehSU0mbNava/Z08bGk/uCG3L+dw42xmNXu94FBCOnThzLbNFOflVrMrKPydUISgjnExOZ/nvzxNfVeTCDjeIQJ7Z8LRD41RN4+uA8tKe2pRKk/tfIqYrBgWdVvEs82fNXtI5/7wAwnPPIva0ZGgjRtw7FUZWlpapOWXFec59cstgtt589i0tngGOFbY46PP8tWUV0m+cpm+Y8cz9K3Z/7XKpxCCSX1DGNu1AV8cj2fh7rg7xMAShn8BXk1hwzPGLOgqeNh6sKznMnI0OUw7Og2DNJjZ7SIi8Hx1PAU7dlKwZ0+1+7fo7oeHvwNHN1xFqymvZu8y8knKdVpOb930H8ehoFDXUYSgDhGfVcwz607hYmfF12Pa4+FgXWmUEnZPN+YJtH0eBi41yxKOy4njiR1PkFmSyZq+axjQoHIfIanXkzZvPmmz38W+cyeCNqzHumFl1FBmYiEb50eRej2fXk81oe+zzbCyMbqX9OXlHP52HZvmz8TW0YnH539Iyz5/PAv4I4QQTHuwKY+3D2D1wRusPnTD/ABrR6N/w9oBvnsMijLMzGHuYUxpO4WjyUdZd7F6dRL3MWOwDmtK2pz30Ofnm9lUahXdR4VSnK8l+tfEaue61vMjtGNXYvbvpqykuJpdQeHvgiIEdYSMAg1Prv0Ng4SvxrTDx/mOmcDuaXDyY2j/Ijy42EwEYjJjeHb3syDgywe+pK1P2wqbobSU5NdeI/frr3F75hn8V61C7VQZP3/1VBo/LTyDlJKhb7SmaadKh3Jeeho/zJrC6a2baNX3AR5fsBQP/8oaRH8WIQTvDWnO4Fa+LNwVx5bzyeYHOPvB6PVQkgM/PmesklqFEaEj6B/Un5XnVnIp65L5tS0t8Z07F31urrE20R34NHSmYbgn5/feprSo+t4GkYMeQVtayoV9u/7ncSko1BUUIagDlGr1PPflabKLtKx9pi2NPB3MDzj4PpxcBe1fggHvVxSMAzibfpaxe8fibOXMVw98RbBrcIWtPCeH2888S+G+X/GePh3vt6ZWJGEZ9AaO/niNvWsv4xXkxPC32+IdVCkQ1347ztdTJ5CbksygSW/R5/lXsLSqMkP5H1GpBIuGt6RdkBtvbrzA6fg7CtDWa2VMNIs/YowmqoIQgpkdZ+Ju4870o9PR6s0f6DZhYbg98zT5m36iNCam2r3bD26IrkzP2V0J1WzeDRoR0Lwl53ZvQxoM1ewKCn8HFCG4z5FS8sbGaC6lFLByVAThd+YJnPwEDr0P4U/AgAVmIvBb6m+8uO9FPG09WTdgHX4OleGf2oQE4keNQnPlCn4rluP25BOVttJytn98geh9ibToWZ/Br4VXlGQwGPQc+e4Ltn44Hze/+jz5wYoaq9ppbaHm06faUN/VlrFfRXEzs8j8gPBR0OZZOLYMYreZmZysnJjdaTY38m+w6vyqatf2eOkl1B4epM9fYO6HANx87Qlp58PFw8loinTVzm3RewCFWZncvnjhrw9SQeE+RBGC+5wVv15ne0wqbw1oQu+m3ubG6B9g11Rj2ehBy81E4FjyMV759RX8HPxYN2AdPvaViVOaK1eIHzUaQ34BAV+sM8sSLszR8NPiMyRdyaXnE03oNiIEtamWf0lBPj8tmM2pLT/Sss8ARsz+oCI3oKZwsbNi3bNtUQnB819GUaC548H8wAfgGwFbXoZ8820zu9bvyiPBj7Du0jpiMs2/+asdHPB6bSKl585RsH1HtftG9AugXGvg4uHqW3E2juyAjb0DFw/u/esDVFC4D1GE4D5mZ0wqS/dd5ZHWfozr1tDceHU3/PwyNOhmDLNUV+YGnk47zcQDEwlyCmJt/7VmtfxLzp0j4amnEVZWBH73HXYRERW2jIQCfnw/isJsDQ+92oqwLpX+gPSb1/l22iSSLsfQ74UJ9B07Hos/2Ju4Jgh0t2f1461JyClh8oZoDFX3UbawhkfXGv0Em180FtKrwhuRb+Bh48F7J99Db9Cb2ZyHDsW6aVMyly1D6swFxt3PgYBm7lw4kES5zvw8CysrmnTpwbVTx9EU3zFLUVD4G6AIwX3KtfRCXt8QTesAF+YPbVG9jPTGZ8GnBYz8zixE9HzGeV759RXqO9Tns36f4WrjWmErPnGC22OeR+3qQtC332DdsEGF7eb5TDYvPovaQsUjU9rg37SyFtGlQ7/yw8wpGAwGRr678J7U7G/f0J3pDzZl7+V0Vh28bm50a2icGcQfMTrIq+Bo5cibbd8kNieWjVc3mtmEWo3nxAnokpLI+/nnaveM6BdAaaGOq7+lV7OFde2JXqfj5tnTf31wCgr3GTUiBEKIAUKIOCHEdSHEW39gXyqEOG96XRVC5FWx6avYttZEf+o6JdpyXvr2LPbWaj55og02llWqaOYnw3cjwNbVGEljXRnLfzn7Mi/vexlPW89qIlC4fz+J417Ays+PoG++MSsXEf1rIjvXxODm58Cjb0Xi7mt0RhsMeg59s5Zdq5ZSL6QJTy5Yhk/jkLv/BzDxbOcgHg73ZcneqxyIMw8bJeIJ45LYr3MgzXwZqH9QfzrU68CKcyvILs02szl0745Nq5ZkrV6N1Jo7lf1CXHDztefSkTuiljDWIHJwc+fab8dqZnAKCvcRf1kIhBBq4GPgASAMGCWECKt6jJRykpQyXEoZDqwEfqpiLv3dJqU032/wH4iUknc2X+RGZhErRkaYZw2XFcH3I4w/R68Hx8p1/+u513lh7ws4WDnwr37/wtOuMpkrf/t2kl6dgHXTpgR+/RUWnp4V9zqx+QZHN16jYStPHn49osIprC0tYcvieUT98hOt+g1k2LQ52Dnf4ai+ywghWPBIS5r4OPHaD+dJziutaoRBK4yCuPkl0OvMznu7/duUlpey9MzSatf0fHUC5Smp5G3aVM0W1sWXjIRCMhMLzW0qFY3bdiT+/Fl0Gk3ND1ZBoRapiRlBO+C6lPKmlFIL/AAM+Q/HjwK+r4H7/i1ZfzqRn84l81rvEDo1rlzbx6CHTWMg/RIMXwc+zStMyUXJjNs7DkuVJf/q9y/qOdSrsBXs2EHKm1Owi4ggYO1a1C4upssZOPD1Fc7uTqBZV1/6j2uOpZVx5lGQmcEPM6dw61wUvZ57kT5jXkJtUTv1CW2t1HzyRGv0BslrP5yjXF/FJ2DvDgOXQHoMnDAvG93QuSFPNn2SrTe2EpdjXnLavnMnbMPDyf58LbLcPCchtL0PagsVl4+mVOtLcLuOlOu0xF84W3MDVFC4D6gJIfADqqZlJpnaqiGECAQaAPurNNsIIaKEECeFEA//u5sIIcaZjovKzKxeG+bvwOWUAmZuvUTXYA/G92psbtw3C67uggcWQnBllE+uJpcX976IRq/h076fEuAUUGEr2LmT5DenYNe6Nf6frkHtYKxOWq7Vs+vTi8QeTyXywSC6jw6tqP6ZcjWWb6e/TkFWJo+8NZuI/g/d/YH/FwLd7Zk3tDmn43NZsf8Of0HTQcbXwfch2zwreUyLMThaOf7hrMBtzHPokpIo3HdHDSN7Sxq19uTqb2mUa82dxn5NwrC0tiEhJrrmBqegcB9wr53FI4EfpZRV/4cFSuOOOaOBZUKIRn90opTyUyllpJQy0tPzP9ewqYtodHpeW38OZ1tLlo4IN99Y5uImOL4SIsdAu7EVzaXlpYzfP57U4lQ+6vURjV0rxaNg126S33gT2/Bw/Nd8gspUSrmsRMfWFee5dSGLriNCaD+4YYUjOvbIATa8+zZWNraMmruYoFat783g/wRDwv0Y1ro+H+2/xsmb5uv+PLgY1NawdYJZFJGztTPjWo7jWMoxTqScMDvFsVcvLAMDjLOCO/IKmnSsh1ajJ+Gi+X3UFpb4N2vB7ZhzNTs4BYVapiaEIBnwr/K+vqntjxjJHctCUspk08+bwEEgovppf38W7orjanoRix5taV5DKO0ibBkP/h2MWcMmyg3lvHnoTS5mXeSDrh/Q2rvyoV2wew/Jkydj26oV/mvWoLI3zgSK88vY/OE50m8V0O+5ZrTsWR8w+Qo2fc+Oj5ZQL6QJo+ctwd2v6kd6fzBnSDMC3e157Yfz5BZXcfQ6+kC/9yDhKJz72uycUU1G4Wvvy9IzS82K0gm1GvdnnkETE0PpmTNm5/iFuGDraMm1qOrRQwHNw8lNTaEgM6OaTUGhrlITQnAaCBZCNBBCWGF82FeL/hFCNAFcgRNV2lyFENam3z2AzsDlGuhTneLotSzWHrvF0x0D6RHqVWkoyYH1j4O1Ezz2JVgYHblSSuaenMuhpENMazeN3oGVG7MU7ttnFIGWLfH/9NOK5aDCHA2bl5wlP7OUga+0JLitMRHMoNez99OVHN/wLWHdevHo9Pfu27167a0tWDkqguziMt7ZYl6JlNZPQWBn2Dfb+HczYaW2YnzEeGJzYtkTb16B1Pnhh1G7uJD9xRdm7Sq1ikatvUiIya5WlTSwRSsAEi6er7FxKSjUNn9ZCKSU5cB4YDcQC2yQUl4SQswRQlSNAhoJ/CDN5+FNgSghRDRwAHhfSvmPEoK8Ei1vbIymkac9bz3QtNJg0MNPY43hoiO+NosQ+iT6EzZd28TYFmMZ0aRyK8aiI0dJmvQ6ts2a4f9ZpQgUZJWyeclZSgu0DJ4QTkCYOwBaTSk/L3qPmP176PDICAa8PAm1xd1JEqspmvs581qfELZfSOWX6CoOXSGM/hNNHhyYb3bOwIYDaeTciDUX1pjNClS2trgMH07RgYPo0tLMzgmO9KZcZyD+QpZZu7t/ILZOziTH/qP+mSr8zakRH4GUcoeUMkRK2UhKOc/UNlNKubXKMbOllG/dcd5xKWULKWUr08/Pa6I/dYkZWy6RVVTG8pER2FpVyRc4+L5xM5YHF4J/u4rmHTd3sCp6FYMbDebViFcr2kvOniNpwgSsGzUyiYAxFyAvvYTNS86iLS1nyKQI6jUybhhfnJfLhnffJv78WfqOHU/nEU/+6dLRtc0L3RrSyt+FGVsuklFYJZTTp7mxBHfU52a5BSqh4oVWL3A97zp7E8zLRLg8Nhz0evJ+NA8lrdfIGXsXa66fMV8CEkJQLziUlGtXan5gCgq1hJJZXIvsiDF+q53UN4Tmfs6VhhsH4PAiCH/cWGTNRHRmNDOOzaC1V2tmd5xd8eDWXLlC4gsvYOnlRcC/PqsoI52TUszmJWfRlxt4+PUIvAJ/b0/i+xlvkJ2cyJA336Fln8q9CeoCFmoVS4a3olSrZ9pPF82dvT2nGXMLdkwxluc20S+wHw2cG1SbFVj5+2PfuTN5P/5oFkoqVIIGrTxIjM2pVnLCN7gJuSlJlBaZ5xooKNRVFCGoJfJKtMzccpEWfs68ULWOUGG6cUnIMxQeXFRRSC61KJWJ+yfiZefFsp7LsFQbl3DKbt3i9pjnUTk4ELBuLRYextyDrKRCNn94FgQ8PKk1HvWNGcjJcbF8P+NNtBoNI2YuoFGbdtRFGns58Gb/UPbFprPpbJXYBFtX6D0Tbh83RluZUKvUjGs5jmu519h/e7/ZtVxGPEZ5WhpFh4+YtQc2d6dcayDlap5Ze73gJgCkXTPPT1BQqKsoQlBLzNl2mbwSHR8Ma4mFqbqn0S/wvDFzePgXYGWK9tEVM37/eLR6LR/3/riidIQuJYXbz40BKQn4/HMsfY1F4jISCoz7CluqGPp6a9x8jde5ceY3fnxvOjYODox+b/E9LRdxN3iucwPaBbnx7i+XyCioskQU8STUC4e9s0BX2T4gaACBToGsubDGbBbh2LMnak8P8tavN7t+/VBXLCxVxMeYh5H6NA5GCJWyPKTwt0ERglrgYFwGP51N5qUejYGI6MIAACAASURBVAjzrRKhc2QJ3DpsnAl4GR3HeoOetw6/xY28GyzuvpiGLsbZQ3l2NrefG4OhsJCAf31WUUAuI6GALcvOY21nwdDJrXHxNuYPXD68ny2L5+EREMio9xbj4lOPuo5KJfjg0ZaUlRt495cqzluVGvrNhYIkOLWmotlCZcGY5mO4knOF39J+q2gXlpa4PDyUoqNHKc+qdA5bWKmp38SVhItZZsJhZWOLm199MhNu3d0BKijcIxQhuMcUlZUzffNFGns5mGcPxx+FgwugxWPGgmomlp9dzsGkg0xtN5VOfp0AMBQXkzjuBXRpafiv+QSbMGNpp8zbhWxdfh4bewsefr01Th62AJzdsYWdH3+If1hzhs+Yh51TFX9EHaeBhz0TewezPSaVfZerxP036ArB/YziWiWcdGDDgbjbuPPFpS/MruM8eBDo9RTs2GnWHtjCg4IsDbmpJWbtHgFBihAo/G1QhOAes3DXFVLyS/lgWEusLUxRQiU5sOl5Y3nlhz6s8Atsv7mddZfWMSJ0BKOajAJA6nQkvTYJTWwsfks/xK5NG8C4wfyWZeewsrFgyKQIHN1skFJybMM3HPjyMxq37cjQqbOxsrWrlXHfTcZ2bUiotyMztlykqKxK3H+f2aApMIqBCSu1FaObjuZY8jGu5V6raLcODsY6rCn5W81TYAKbG0Ntb182Xx7yDAiiIDND2dRe4W+BIgT3kHO3c/n6ZAJPdwyiTaCpRLSU8MsEKM4ybrhiKit9JecKs4/Ppo13G6a2m2o6VJI6cxbFR47gM3sWjj17ApCVVMTWZeextFbz8OsROLnbIg0G9q/7hJObfqB5z34MmvQWFlZWtTLuu42VhYoFw1qQVqBh8e4qDlzvZsbIq1OfQm7lfsSPhTyGrYUtX13+yuw6zoMGo7l4kbKbNyvaHN1scPa0JfkOh7FnoHEpLvN2fM0PSEHhHqMIwT1Cb5DM2HIRL0dr3ugfWmk4/x3E/gK9Zxg3aAfyNHm8duA1nKydWNx9MZYqY4RQ1sqV5G/ejMcrr+D62GMAZCcXsWXZOdSWKqMIeNiiLy9nx0dLOL97O5GDHqHfC6+iUqur9envROsAV57qEMiXJ+I5n1jlod1zGggVHJhX0eRi48KQRkPYdnMbmSWVBQydBj4IKlW1WYFfE1dSruWZ7ZT2uxBkJcTflfEoKNxLFCG4R3z7WwIXkwuY8VAYDtamks45t2DnFAjsAh3HA0bn8NQjU8koyWBpj6UV20zm/rCerFWrcX50GB7jXzGenlrMlmXnUKkFD0+KwNnTDl2Zhi2L53Ll2CG6jn6G7k88V2cSxf4qb/QPxdvRhrd/iqksV+3sB+1fhAsbIL3SofxU2FMYpIHvr1SWvrL08sK+Y0cKftlm5hyuH+KKtrScrCp7FDi4uWNj70BWYvxdH5eCwt1GEYJ7QEahhkW74+jS2IOBLUzROga9cc9doYahnxgjXYCV51ZyPOU409tPp6VnSwAK9x8gbc4c7Lt3o96sWQghyE0rZsvScwhhFAEXbzvKSorZNH8mt86foe+48bQb8mhtDblWcLSxZNagMGJTC/ju1O1KQ+eJYOUAhyqL9vk7+dOtfjc2XduEVl9ZwM7pgQHokpPRXKoUDd8Q4x4OSVdyK9qEELj6+pGb+u/qKyoo1B0UIbgHLNhxhTKdgTlDmlV+Oz+6FBJPwsDF4GKs9Lknfg+fX/ycR0MeZVjIMABKz58n+fXXsQkLo/6HHyIsLcnPLGXL0nNIKRkyKQJXH3tKiwrZ+N47pF6L46GJU2nZu25lC9cUA5r70KWxB4t3x5FdVGZstHODji/D5S1mpSdGho4kR5PDvoTKPQkcevcGtZrCPZUF6uydrXGtZ09yXKUQALj6+JKbmnp3B6SgcA9QhOAuc/JmNpvPJTOuW0Maehrr/5B81hgq2nwYtBgOGLeafOfYO7T0bMnb7d4GQHv7NokvvYyFl5dxTwF7e4pyy9iy7Bzl5QaGvBaBWz17SvLz2Pju22QlxjPkjXcI7diltoZb6wghmD04jBKtnoW7qjiOO7wM1s7GGk4mOvp2JMAxgPVxlYlkFq6u2LVrS+GePWbLQ37BLqTezDfzE7jW86MwOxNdmbJ1pULdRhGCu4hOb2DGzxep72rLKz1NOQO6Utj8Ajh4G7dZFIISXQmTDk7CzsKOD7t/iJXaCn1BAYkvvgQGA/5rPsHC3Z3SQi1bl59DU6xj0KvhuPs5UJSTzfp33yY3LZWhU2bRsHXb2h30fUBjL0ee69KA9VGJlY5jWxfo+Apc2QYpxhLSKqHisdDHOJtx1mw7S6d+/dDGx1N2rTK81KehEzqNntzUynBRl3rGTO68dPPKpQoKdQ1FCO4iXx6P51pGEbMHNausLHpgPmRdhSEfg60rUkrePfEutwtvs7DbQrztvY25AhMnok1MxG/lCqwbNKCstJytK85TkK3hoVda4h3kREFWButnv0VhdhbD3n6XwJbhtTvg+4gJvYPxcrRm5paLld/iO7wINi5ms4KHGz+MtdqaDXEbKtoc+/QBISjcU1mp1LuhMQkv7WZ+RZurj0kIUqvvb6ygUJdQhOAukV1UxvJfr9E9xJM+YcZNYEg8bdxkvc2z0MiYA7Dx6kZ23NrBK+Gv0K5eO6SUpM15j5ITJ6n37rvYt2uHTqtn+8fR5KQU88ALLfANdiUvLZX1s9+itLCAR6e/R/2w5v+hN/88HKwtmD6wKReS8lkfZdpS28YZOr0KV3dCsnFXMmdrZx5o8AC/3PyFIm0RABaenti2bm3mJ3D2tMXGwdJcCEwzgtw0RQgU6jaKENwllu67SolWzzsDTZvN6DSw5WVw8oO+cwCIzY7lg1Mf0Nm3M8+3eB6AnHVfkLdxI+7jxuHyyFD0OgM7P4kh7UY+fZ9rRmBzd7KTE1k/eypajYbhM+bhG9KktoZ5XzO4lS/tGrixcNcV8kt1xsb2LxgrlB6uzDYeGTqS0vJStt/cXtHm1K8vZVevok0wJqIJIfBp6EzazYKKY6xs7bBxcKQgqzIXQUGhLlIjQiCEGCCEiBNCXBdCvPUH9meEEJlCiPOm1/NVbE8LIa6ZXk/XRH9qm7i0Qr777TZPdggk2NuYKczBBcYlocErwMaJQm0hkw9NxsXGhfld56MSKgp//ZWMRYtw7N8fz9cmYtAb2LP2EomXc+jxRBMat/Ei83Y8G959G4PBwIiZ8/Fu2Pg/d+YfjBCCWYPCyCvV8fGB68ZGa0djXkHc9oq8gjD3MEJdQ9l8fXPFuQ6mrO2iQ4cq2nwaOpGXXoKmWFfR5ujhSWGWsn+xQt3mLwuBEEINfAw8AIQBo4QQYX9w6HopZbjp9S/TuW7ALKA90A6YJYRw/at9qk2klLy37TKONpZM7B1sbEw6A8dXQOunoVEvpJTMOj6LlKIUFndfjJuNG6WXLpH8xpvYNG+O7/sLAMGBr69w81wmXYYHE9bZl/Sb19nw7tuo1GpGzH4fj4Cg2hxqnaCZrzPD29Tni2PxJGSbHL3txoGlvTGEF6NgDA0eyqXsSxVOY6uAAKwaNKDo0OGKa/k0qO4ncHT3oFCZESjUcWpiRtAOuC6lvCml1AI/AEP+5Ln9gb1SyhwpZS6wF6jTAfC/xmZw9HoWk/oE42pvZVwS+vklcPQ1lkYGvrvyHXsT9jKx9UQivCLQpaeT9NLLqF1d8F/1McLGhiMbr3HlZBrtBjWgVW9/Uq/HsfG96VjZ2jFi9ge4+dav5ZHWHSb3C8VCLfhgl2n/ADs3iHzWuHFNjrGC6MAGA7FUWfLz9Z8rznPo3p2SU6cwFBsFxDPQOLurmmHs6O5JYbb5vsYKCnWNmhACPyCxyvskU9udDBNCXBBC/CiE8P8fz0UIMU4IESWEiMrMvD+/gWnLDczbEUtjLwce7xBobDz0PmTFweDlYONETGYMi6MW06N+D55u9jSGkhKSXnoZQ1ER/qtXY+HpSdSOeGIOJBHex5/IB4NIvR7HpnkzsXF0ZMTsBbh4+/znjiiY4e1kw4vdG7EjJo3T8aaS1B3HG7O5j68AjPWHegX0YtvNbej0xqUfhx7dkTodxSdPAmBlY4GLtx2Zt4sqru3k4YmmuAhtqXmZagWFusS9chb/AgRJKVti/Nb/5f96ASnlp1LKSCllpKenZ413sCb46kQ8t7KKeWdgUyzVKmNkyrHlxh2zGvehUFvIm4ffxMvWi7ld5iIQpEybXlFS2iY0lIuHkzn1yy2adPCh07DGpN24WiECj82cj5OHV20Ps04ytmtDfJxsmLvtsjGc1KkehI+Gc99CoTEPYGjjoeSV5XEw6SAAdq1bG5P4Dlb6CTz9Hci8XXVGYKwFpcwKFOoyNSEEyYB/lff1TW0VSCmzpZSmfH/+BbT5s+fWFfJKtKwwhYv2CPUCvQ62TjAmjvWfh5SSOSfmkFacxgfdPsDZ2pnsNWso3LULrzcm49C9OzfOZnDo+zgCW7jT48kmigjUILZWaqYMCCU6KZ+t0aZwz84TwaCDEx8D0KFeB7ztvPnp2k8ACCsr7Dt3pujw4YosY48ARwpzNGiKjLMGRw/jlxLFT6BQl6kJITgNBAshGgghrICRgFkdXyFE1X0RBwOxpt93A/2EEK4mJ3E/U1ud4+MD1ykqK2fag6Zw0RMfQfpFeHAx2Djz8/Wf2RW/i1fCXyHcK5zC/fvJXLYcp8GDcHvuOZLjctmz9hI+DZzpP7Y5GbeuVRGBBYoI1AAPh/vRws+ZD3ZdoVSrN24E1OwRiFoLpbmoVWqGNB7C8ZTjpBUbZwkO3btTnp5O2RWjf8EzwOgnyDT5CRzdjDOCAmVGoFCH+ctCIKUsB8ZjfIDHAhuklJeEEHOEEINNh00QQlwSQkQDE4BnTOfmAO9hFJPTwBxTW50iKbeEL48nMKx1fUJ9HCHnpjF7tclD0PQhbuXfYsGpBbT1actzzZ+j7No1UkwRQvXmzCErqYgdqy/g7GnHwFdaknX7Oj/OnVFFBO7PpbC6hkoleGdgU1LzNXx+1LT5TJfXQFsEZ74A4OFGD2OQhoqcAoduXQEoOnoUAE9/kxCYlofsXIyVSUsLKiOJFBTqGjXiI5BS7pBShkgpG0kp55naZkopt5p+f1tK2UxK2UpK2VNKeaXKuWullI1Nr3U10Z97zYd7riIETOobYtxxbNskUFnCg4vQ6rVMPTwVK7UVC7osQOYXkPjyKwh7O+p//BGFhQZ+WRmNla0Fgye0Ijf1Jj/OnYGtk5MiAneB9g3d6RfmzZpDN8kp1oJPC2jQHX5bA+Va/J38aeXZiu23jEJg4emJdUgIJSdOAGBjb4mjm03FjMDSyhorW1uK83P/7T0VFO53lMziv8jllAI2n0/mmc5B+LrYwoX1cPMg9JkFTr4sP7uc2JxY5nSag5e1O8mvv055Whr+K1eis3Vl64poDHoDgyaEU5idoIjAPeDN/qEUa8tZ9XuSWadXoTAVLhkTygY2HMi13GsVOQX2HTtSEnUGg8ZYZdTD34GsxMrIITtnF0ryzLeyVFCoSyhC8Bf5YNcVnGwsebl7YyjOhl1vQ/12EDmGo8lH+eryV4wMHUmvgF6kf7CQkhMn8Xn3XdShzfll5XlK8st46JVWlBUn8ePcGdg5OSsicJcJ9nZkWOv6fHUigeS8UmjUGzxC4cRKkJL+Qf1RC3XFrMC+U0ekVkvp2bMAuPs5kJ9RQrlOD4CdsyslytKQQh1GEYK/wPHrWRy6msn4no1xtrOEPdOhrAAGLSerLIfpR6fT2KUxkyMnk7dpE7lff43b00/hOGgIOz6JISe5mAHjWiANaRUiMHzmfEUE7gGT+oaAgKV7r4JKZSxRnRYD8Udws3Gjs19ndt7aiUEasIuMBEtLik3LQ26+9kgJeenG3AE7J2dK8pUZgULdRRGC/ycGg2TBziv4udjyZMdAuHEAor+Hzq9h8GrCO0ffoVhXzKJuizBciCV19rvYd+qEx+Q32LvuEslxufR6qgnWNtmKCNQCvi62PN0xkJ/OJnE1vRBajgA7Dzj+EWDMNE4rTuNM+hlU9vbYtWpF8bHjgFEIALKTjRnH9i4uFCtCoFCHUYTg/8n2mFRikvN5vW8INrIMtr0Gbo2g25t8fflrjqUcY0rbKQSWOZA0YQKW9erhu2Qxxzbf4sbZTDoNa4yLdwk/zjf6BBQRuPe83KMx9lYWxp3MLG2g3Vi4thsyr9LDvwe2FrYV0UN2nTqiiY2lPDcXFy87VCpBjmmTGjtnFzSFBejLy2tzOAoK/28UIfh/oC03sGh3HE18HHk4wg8OL4TceBi0jNiCWyw7u4xe/r0YFjSEpIkTkSUl+K/6mJioQmIOJNGqtz/1Qw38OG8GVrZ2PDZDEYHawNXeihd7NGJfbDpR8TkQOQbU1nByFXaWdvQJ6MOehD1o9VocOnUCKSk5eRK1hQoXHztyUn4XAmOdxNLCgv90OwWF+xZFCP4frI9K5HZOCVMHNEGdFQfHV0Kr0Wj82/HWkbdwtXbl3U7vkjFvPproC9R7fwG38xw5vuk6jVp70bSDNT/OfQcLS0ujCHgqyWK1xbOdg/B0tOaDXVeQ9h7QaqRxia84i4ENB1KoLeRI0hFsmjdH5eBA8fHK5aGcFGPkkI2DcS9qTVHhv72PgsL9jCIE/yManZ6P9l+jbZArPUI8YMcbYGUPfeew/OxybubfZG7nuchf9pK3YQPuY8dS1KAt+764TL3GzrR5wJmN894BYPiMebj41Psvd1S4m9hZWTChdzCn43M5EJdh3OS+XANn1tG+XnvcbNzYfms7wsICu3btKD51CgC3evYUZGnQlemxcTAmmWmKi/7TrRQU7lsUIfgf+eZkAukFZUzuF4q4uAnij0DvWRwvuMY3sd8wusloIrIdSJ/zHvadOqEe8Tw7Vsfg7GFLl0e9+WnBDAzl5QyfMU8pJX2fMLKtP0HudizcFYfBIxQa9oDTa7GQkr6BfTmSdIQSXQl2bduiS7iNLj0Dd1/jLCA3rRgb+99nBIoQKNRNFCH4HyguK2fVwRt0aexBB18L2D0NfCPIbz6UGUdn0MC5Aa8GPUXShIlYeHriPGsB21bFoLJQ0WO0L1sWz6Jco+HRd+bi4R9Y28NRMGGpVjGpbwhX0grZcTEV2r0AhSlwZRv9g/qj0Ws4nHzYGEYKlESdNosc+l0IypQZgUIdRRGC/4EvjseTU6xlcr8QODAfijORDy7hvVPzydHk8H7HeWRPmY4+NxevJcvY9U08muJyej3hz86P3kNTVMij78zFK6hhbQ9F4Q4eaulLsJcDS/deRd+4H7gEwKnPaO3VGncbd/bE78GmaRNU9vaUnD6Nk6ctKgtBXnox1g7KjEChbqMIwZ8kv1THmkM36NPUiwjL23DqU4h8jm1lKeyO383L4S/j/uUuSn77Da+Zszl0RE92cjE9RvlzYN0CivPzGDZtjrLH8H2KWiV4vW8INzKL2XIhDdqOhYRjqDNiK5aHSqUW2zatKTkdhUolcPawJS+9FGs7O0DxESjUXRQh+JP868hNCjTlTOrTGLZPBls3UjuMY/5v8wn3DGd4oi85a9fiMmo054tDuX05h05DfTn54xIKsjN55K1Z+IY0qe1hKPwH+jfzIayeE8t/vYau1eNgYQun1lQuDyUdxi6yLdobNyjPzsbF2468jBJUKjXWdvbK0pBCnUURgj9BdlEZa4/eYmCLejRL2wpJpzH0ncP0MwsxSAPv1RtL2oyZ2EZEkBg+itjjqYT38eLCvo/JS0tl6JSZ1G/avLaHofBfUJlmBQnZJfwUWwwth8OFjUQ4BuJh68Hu+N3Ytf3dT3AGFy878jNKMRgk1vYOyoxAoc6iCMGfYM3hm5Tq9Ezu4g77ZkFAJ76y1HI67TTTwiaif2s+Kns7ip6exemdtwmOdObW2c/JSbrNkDemE9C8VW0PQeFP0rupF638XVjx63V0bcZCeSnq898Zl4eSjyBDGyJsbSk5fRoXbzv05QaKcjTY2DsoMwKFOkuNCIEQYoAQIk4IcV0I8dYf2F8XQlw2bV7/qxAisIpNL4Q4b3ptvfPc2iajQMOXx+N5OMKPhtGLQVNAXLcJrDi3kl5+PWn96VG0ycmIyYs4vDUFvxA7shK+JzP+JoNef5ug8Db//SYK9w1CGGcFyXml/JDoDIGd4fS/6B/QlzJ9GYfTj2MXEV4hBGAsPmfjYE+pklCmUEf5y0IghFADHwMPAGHAKCFE2B2HnQMiTZvX/wgsrGIrlVKGm16Duc/4+MB19AbJm80K4OxXlLUfx9uXP8PJyonJlwMpOnAAy/HvcOBAGS4+VmjyNpN2/SoPTZxKozbta7v7Cv8PugV7EBnoysf7r6Nt8zzk3SYiLw1PW092x+/GNjKSsqtXcbQx7lucm16Ctb0DZUrUkEIdpSZmBO2A61LKm1JKLfADMKTqAVLKA1LKEtPbkxg3qb/vScot4btTtxnRph71jkwDR19WujhxLfca71s+RvHqz7Ec9BhHbvlhZQsW7CL56mUeGD+Z4Padarv7Cv9PhBC83i+EtAIN3+a1ACc/VKf/Rd/AvhxNPooqogVIiYy7gJWNmvz0Eqzt7NFqSmu76woK/y9qQgj8gMQq75NMbf+OMcDOKu9thBBRQoiTQoiH/91JQohxpuOiMjMz/1qP/yQf7b+OQPCmx3FIiyGq8zi+ivuBp90exO2Dr1A3aU6U4wPoNDocHA6ReOkc/ca9StPO3e9J/xTuHp0aedCpkTsfH45HG/E03DxIf5cwyvRlnHLNQVhaUnrubEXkkKWNDdpSRQgU6ib31FkshHgCiAQWVWkOlFJGAqOBZUKIRn90rpTyUyllpJQy0tPz7lfqTMwp4cczSYxp7YjLiYWUNOjCO8m7CbDxZeiXN9DrJVc6vU52aglu3idJiDlFj6fG0qJXv7veN4V7w+R+IWQVafm+vCeoLAm/eQJ3G3f2pR3GplkzSs+ew9nLjtz0Eqxs7NBpNEgpa7vbCgr/MzUhBMmAf5X39U1tZggh+gDTgcFSyrLf26WUyaafN4GDQEQN9Okvs+rgdVRCMF6sh7JCPgxoQkpRCvPPNabs4iWSHpvH7evFePmdISH6GJ0ee5w2A4f89wsr1BnaBLrRPcSTZSfy0IU+hCr6O3rW78bR5KNYhbdEc/EiLh7WFOWUoba0QkoD5dqy/35hBYX7jJoQgtNAsBCigRDCChgJmEX/CCEigDUYRSCjSrurEMLa9LsH0Bm4XAN9+ksk5ZawMSqJ15qVYn/ha06EP8L623uYmtcRy59/JeeRqcTdELj7XOB2zCEiBz1Ch0dG1na3Fe4Cr/UJJrdExw6rB0CTTx+9JSXlJSQE2iJ1Ouy02QDodRYA6Ewb3Cso1CX+shBIKcuB8cBuIBbYIKW8JISYI4T4PQpoEeAAbLwjTLQpECWEiAYOAO9LKWtdCFYdvIEQkjFFn1Bo78ZM7S3alfrS5ovT5LcfSnRuAE5ul0i58iut+j5It8efRQhR291WuAtEBLjSLcSTd2NcMbgH0y7uAI6Wjux3Nk56rVOvAaDTGv8rKX4ChbpIjfgIpJQ7pJQhUspGUsp5praZUsqtpt/7SCm97wwTlVIel1K2kFK2Mv38vCb681dIzitlY1QicxvFYZ3yGwtD2pOfn8nkzXqK3BoR7dQXW9vLZN7cQ1jXnvR+7kVFBP7mTOwdTE6JjhOuQ7BMPkNX9xbsKjiJpb8/6qtnAdBqTEKgRA4p1EGUzOI7WH3wOnZoeDR7DYf8wvg5N4b3jwVSllHKhZYvg7xCbvIuGrftSP+XXkOolD/h3502ga50DfZg+s3mSAtb+hQVkFeWR0nTAMrPncLa1oKyUuOXAWVpSKEuojzFqpCaX8qG00ks891PUWkGsx3UPHHFA89TSVzq/g6akqsUZ+0gqFVrBk6cgkqtru0uK9wjJvYOJr7Eique/egcdwhrtRWX6pWjz87G0VmNxpRUrMwIFOoiihBU4ZODN/CXqfTIXs/8RuG4JxTy0LZcrnR5g9yiZDQFO/ANDWPw5GlYWFrWdncV7iGRQW50buzO/IyO2GmL6Wjrx07HeADsKKKk0Bg2qlOEQKEOogiBifQCDd+fTmSl+4/sc7DncEEm07ZZc6P5E6SW6ygv2YZXUEOGTp2FpbVNbXdXoRaY2DuEQ8X+ZDk2pU92Kufts5D2dtgUpFKcbxQCxVmsUBdRhMDE6oM36CLP4V18grkeHv/H3nmHR1V8f/id3Ww22fTeGySE3nvvCEi3YPnaERUQFAURQRFBsFBUUECKCla6Si/SpAUIEFoCKZCQ3nvZnd8fuyRZwPYjIZT7Ps8+e3dm7r0zeeB+7syZcw6TtupItWlDnK0X+qKNOHl5M+ztaeVJSBTuP1oHOdOuliuLCrrQNekiKrWa9NouaBIikQbjDLFEsREo3IUoQoAxwujqI9HMsvmeD7z96HGgGMe0ICL92lNWuB57VxceeucDrO3sa7qrCjXM2J4hrMpvjbWwpqXKluPu+WgunwVhFAJlaUjhbkQRAuCrPdE8wWaOatJJjtfT+7g3EfWGUla4Dp2DLQ+/MwNbJ+ea7qbCHUDbWi40CvJmg+xEj9R4DrllYV2YhhAWCJVaMRYr3JXc90KQklvEtsPhPGm9noVWboza6sKpxk9RnL8erbWGh9+Zgb2be013U+EOYmzPEJYUdqN7Xi4XvQSWxVkIJGoLrWIjULgrue+FYPGeaMaK75ntYMMLGzVcCP4feYWbsNAaeHjKBzh7/10gVYX7kXa1XHAMaEKioQ7BakGKtyXWhlyESqP4ESjcldzXQpCaW8zpwzuwcDhOwGEt+U5PkKXfh1pdzEOT38ctIKimu6hwByKEYGzPEJYXd6NHThanPIrQ5iYipYbSYkUIFO4+7mshWLL3Ii9qlrEj1YlaWUNIsYhAJ4+bkQAAIABJREFUkMvQSe/iFRxa091TuINpX9uFZN8+tC4QXPAR6PJSkAa1IgQKdyX3rRCk5RWTf3gFWzT59AnvTJxdMhjSGfzmZPzqN6rp7inc4QghGNWrAQcLO1LooceqKA2D1FBcoNgIFO4+7lsh+Hb3SWrbrKX1vuZcchVIfSL9x75JULOWNd01hbuEjsGunHIfTDOLAgwyDYGG4vyCfz5RQeEO474Ugoz8EtQnPqT4eDAJTu4Yyi7T68VXCW3Xsaa7pnAXIYTgkQe64ZLnRaJ9GggNxQXK0pDC3cd9KQSrt2yjLCmJUnVd9GUxdHp8BI179KrpbinchXQOcSVeN4RktzSE0CgOZQp3JfedEGTmFVN4ZAbazNaU6i/RvN9wWg9SUkwq/P8QQtD2gf9h41yAyiApK1aEQOHuo0qEQAjxgBDighDiohDirZvUa4UQP5nqDwshAivVTTKVXxBC9KmK/vwda76fhUVcC4oNl6jTujfdnn6yum+pcI/TpZ43xc6NsdCXodeX1nR3FBT+M7csBEIINbAA6AvUBx4TQtS/rtnzQKaUMhiYC8w2nVsfY47jBsADwELT9aqFlPRUSnYnUCRj8QxoyoOvj6muWyncRwghaNv1TaAYMGDQ62u6SwoK/4mqmBG0Bi5KKaOllCXAj8D1ay2DgG9Mx6uBHsKY33EQ8KOUslhKGQNcNF2vWlj9+iQKDVextvbi8VnvKykmFaqMTi1bIdXFAGSlJNVwbxTuRWLOR/DV8yPJzcmu8mtXhRD4AFcq/Y43ld20jSnZfTbg8i/PBUAI8aIQIkwIEZaamvqfO6nX69GXSjQqN15cslBJMalQpQgh0DgaQ5RvWj6nhnujcK9RXFTMxmkfk5+XwJYvP6/y61tU+RWrCSnlYmAxQMuWLeV/PV+tVjNq5RLyszKU7GIK1ULtZq0J3x5N3sVY9AaJWqXMOBVuHX1ZGV+/PJoyQzrWWjcenvhOld+jKl6LEwC/Sr99TWU3bSOEsAAcgPR/eW6VoVKpsHN2ra7LK9znuNfyB0CVb8uOAwdruDcK9wIGg56Vb75NUUEiOuFLh7der5b7VIUQHAVChBBBQghLjMbfjde12Qg8bTp+CNglpZSm8uGmXUVBQAhwpAr6pKBw27F1tgXAusSeM4fmoTf854mrgkI5UkrWzphN2tWz2BCMoXkZTaop/M0tC4FpzX80sBU4B/wspTwjhHhfCDHQ1Gwp4CKEuAi8DrxlOvcM8DNwFtgCjJJSKlsuFO5KtDprAKTKFnXxcTadiK3ZDinctUgp+f3zL4iL+BNrURdD4G7aPTSj2u5XJTYCKeUmYNN1ZVMrHRcBD//FuTOA6huhgsJtQqO1AqDY0oa0LDV5O76jf7MpqBRbgcJ/ZNeK5Vw4sBVLdX0K7XdT4v8CzWt5VNv9lK0zCgpVRIUQWOOQpKZByW9sikis4V4p3G0c+PkHwresxUJTH21pGG5eVvTq/0S13lMRAgWFKkJjZRQCvUpFcKIV6bZXWbdtFwbFVqDwLwn7bT2H1qxCramLf8oF1K0ucdzjdVoGVm/OdEUIFBSqiGtCgCzBI9eVvZY2dMj6la1nFAczhX/m1M4t7Pnua1SaYEIT09nT9Ry5OV14vE+Xar+3IgQKClWExlILgJSllGhdyc/U0FW7ny93RCizAoW/5dz+P9i+eAEqi0CC0y34veUR/mfQE+77HK2Dqnc2AIoQKChUGUKlwsLSEmQpBdauhCZIwnQQkrqDbWeTa7p7CncoUUf+ZPMXc1BZ+OBT4E2U8ybaeWSxvOAxXup9e7IlKkKgoFCFaLRWqDV6SlyDaJpkxU5HF56z2s1nO6Mwus4oKFQQE36M3+bNRqg9cJFNsM79mbMdJXVy/Un170fbWi63pR+KECgoVCEaK2s0lgaKnXwJSZAcUYMfkcik0+w4l1LT3VO4g7hy9jQbPvkAoXLB3qoboZe+YeEgydT0VKYUPcXYXnVuW18UIVBQqEI0Wi1qCwMFGke0OYW4ZBnYa2PPSJs9zN8ZqcwKFABIjLrAulnTQDhgZTOAJie+4rMHcnm2JJ3Dhp7YBzSl3W2aDYAiBAoKVYrGygqVuoyCYjUGoaJFqi27PGvRX+4lJiGZ3ReUWcH9TkpsNGtmTkVijUY3hCbhS9nWMhvbIB0DitR8UDCEsT1DbmuYfEUIFBSqEI3WCkEpUkKJkw8dM9zYL/PQGwp52j6M+TsUW8H9THrCFVbPmILBYIFKO5QGkevI8Erjt44WTLtyic/ko4QE+NG+9u2bDYAiBAoKVYpGq8UYfgsM9VoReLmYQn0xhz3r8ILVH5yMz+KPyP+eT0Ph7icrOYnV0ydTVmJAWA6lTuYJHEvPMaNPHhNyi7DUhbAorxOv9ri9swFQhEBBoUrRWFkjDSUAlPrXwyImAVeDDTs9gnDOOUcP+wRlVnAfkpuexi/TJ1NcWIywHIy/yMLvwjpmDiihubMPg1LjmVz8P5r4O9Mp5PaHyleEQEGhCtForSgrLUZlISh28gWDgYEl9fij6Cp6jQ2T3A8SfiWLfVFpNd1VhdtEflYmv0yfTEFONirtYDztdQTvm8vmQV4k+Fnx7sWTxHr3Z3NOEGNrYDYAihAoKFQpGistpUVF2LtYU6B2ACFom+ZEZnEWJ+r2pHbSFurY65mv+BXcFxTm5bJmxhRy0lLR2g3F0dGD0G3vkdY5lBXBV5lscMIVFePSh9DEz5EuddxqpJ+KECgoVCEarRWlxUXYuViRk1WKNjgYz5hsLFWW7HR2R5QV8kHtsxyLy+TAxfSa7q5CNVJcUMDaD98l/Wo8dm7D0Gh9aHDwYzQB7kxsE00v1yb0vXiAU7Ve4GS2jnE1NBsARQgUFKoUSytrpMGAvbOGnNRCrJo0peRUBG0927A7IwLp3ZSWaevwtNMqfgX3MKXFRaz/6H2Soy/iGvAQpWU+NEv4GV1ROvOGadDq7Jgcex6ca/Ha5Q409nWga2jNzAbgFoVACOEshNguhIgyfTvdpE1TIcRBIcQZIcQpIcSjlepWCCFihBDhpk/TW+mPgkJNo9EaA8/ZOKooKdKjbtQcQ04OfdWNSchL4EKDB1Glnufdpjkcjc3kYLQyK7jXKCstZcMnM4g/fwafBo+Sk+FNS+vTWJ/cxfEXO/CnKob3nFvhknaRvbXGE51Zxrjb7DdwPbc6I3gL2CmlDAF2mn5fTwHwlJSyAfAAME8I4Vip/k0pZVPTJ/wW+6OgUKNYmJLTWNsZ/1OX+tYFoGmSNSqhYqeVBrT29C7chIe9lvk7omqsrwpVj76sjN/mzSbu1AlqtxpOWrwXzUIKsdm4kLInBzFLu4shAX3oFvYj+uBeTDrtRVM/R7qFutdov29VCAYB35iOvwEGX99AShkppYwyHV8FUoCamwMpKFQj13ISXBOCAo0jagcHVGciaerWlF0J+6Dxo6jPbWRsO2cOx2RwSJkV3BMYDHq2LJzLpbBD1OsynIQoL0IbWuP07WSs2rRmYmg4XjZeTEhJAX0pG7zGcjW7iPG969TobABuXQg8pJTXcvElAX+bVFMI0RqwBC5VKp5hWjKaK4TQ/s25LwohwoQQYampikOOwp2JpUkIrHQSBGSlFGLdrBmF4Sfp7t+dyMxIrtTvB/piHrbYh5udlnk7Imu41wq3ipSSHUsWcP7AHhr1eJi4CB98Q+zxX/8eFk5O/PS4N/EFV5lR6yFsz26gtP04Zh0qpnWgMx2Db7/fwPX8oxAIIXYIISJu8hlUuZ00Wr3+0vIlhPACvgOelVIaTMWTgLpAK8AZmPhX50spF0spW0opW7q5KRMKhTsTjdYaAH1ZMfYuVmQl5WPdtCklly7R3aEVALsKroBfWzQnvuGVLkEcis7gwEXFr+BuRUrJH98s4fSubTTqOZjYs4E4eupoGLkcfXIiqZOfYVXSbzxT/ylaHFgEToGsVA8mJbeY1++A2QD8CyGQUvaUUja8yWcDkGx6wF970N80opYQwh74HZgspTxU6dqJ0kgxsBxoXRWDUlCoKaxsbQEozs/H0UNHZnIB1s2aAeAUnUaoUyg7L++Els9CxiWecI/F28GKj7deUHYQ3aUc+GklxzdvpGG3/iREhWJpZUF7+5OU7NmF/Ruv8nbGcuo41WF0EZAWSVGvWSzYn0CHYJfblm/gn7jVpaGNwNOm46eBDdc3EEJYAuuAb6WUq6+ruyYiAqN9IeIW+6OgUKNY2RiFoCg/DycPG7KSC7Bq0ADUagqOH6dnQE9OpJwgKbAd6FywPLaUMT1CCL+Sxa7zSmTSu41Da37k8LqfqN+pJ+nJzSgrMdCjg4GCRfOwH/Agn/hFkFOSw8wmY7Hc+ynUfZDlKXVIyyvh9V6hNd39cm5VCGYBvYQQUUBP02+EEC2FEF+b2jwCdAaeuck20VVCiNPAacAV+OAW+6OgUKNoKwmBo6eOshIDBSVqrOrXpyAsjAcCHwBgW/weaPEsXNjEQ7XKCHDR8cm2SCW38V3EkQ2rOfDzSup17EZRSUeyUwrp9bA3hTMmoA0O5vgzbdhxZSejm40m9NASkJK8btNZtPcSXUPdaBFww277GuOWhEBKmS6l7CGlDDEtIWWYysOklC+YjldKKTWVtoiWbxOVUnaXUjYyLTU9KaXMu/UhKSjUHJbW1giViuL8PBw9dABkJRWga92KopOn8Lf0oJ5zPbbGboVWz4NQoTm2lHE9QziXmMPmiKQaHoHCv+HY7xvY9/0KQtt1wsK6F1cjc+j2WDByziRkWRmWs6cw8+SnNHdvztNaXzi3ETq/wdLTerIKSnn9NmYf+zconsUKClWIEAKtjS1FeXk4eRqFIDO5AJs2bZClpRSGh9MnsA+n0k4RLwxQfxAc/46B9RwIcbdlzvYL6JVZwR1N+Nbf+ePbJYS0bo+T3xCiwtJoMzAI+18XUHTmDF4fzeLdK19hkAZmtHsX9eaJ4BJMdtOX+Hp/NL3re9DY1/Gfb3QbUYRAQaGKsbKxoSgvF529JRorNVnJBVg3bwFqNfmHD9MnsA+AcVbQ9mUozkZ96kde71WHS6n5rD+RUMMjUPgrTu3cys5lX1KrRWsCmz9B+PZ4GnTyJij9T7I3bMB19GjWeyRwJOkIE1tPxPfkGsiIhn4fs+RgArlFZbx2h80GQBECBYUqx8rWjqL8PIQQOHnoyEzKR21rg1XDBhQcOYqvnS+NXRsbhcC3FXg3gyOLeaCBOw197Jm3M5KSMsM/30jhtnJmz062L/mCwKYtaNRzBPt/vkRgIxdaBueQMnsWtt27k/lYT+Yfn09X364McWkG+z6B+oPJ8OzI8gMx9G/sRT0v+5oeyg0oQqCgUMVY2dhSnG80dzl66shKLgDApnVrCk+fxlBQQJ/APpzLOEdc7mVo8zKkRSKidzO+dyhXMgr55diVmhyCwnWcO7CHrV/Ox79BY9oOG8PO5Rdw87ej24OuXH3tNSz9/HCZ8R4T97+FnaUd77V/D7FlEgg19JnJF7suUliq57WeITU9lJuiCIGCQhWjtbGlyCQETh425GUWU1JUhq51GygtpeDECXoH9gZgS8wWaDAYbNzh8CK61jHuJvl850WKSvU1OQwFE5GHD7D5i0/xqVuf7s+9wdYl59A5WNL3hbokv/EasrgY3wVf8FnUUi5mXeSDjh/gEvsnRG6GrhOJNzix8lAcD7fwI9jdrqaHc1MUIVBQqGKsbGwpys8HwMnLZDBOKkDXvBlYWFBw+AieNp40d2/OltgtYKGFls9B1DZERjRv9A4lKaeIVYcv1+QwFICLYYf5ff5HeAWH8sCoSWxedB6DQdJ/VBNy5s6i6PRpvGfP4ohlAqvOreLJek/S0aUxbJoAHo2g7SvM3R4FAsb1ujNnA6AIgYJClWNla0dxXh7SYMDFx+hXkJ6Qh8rGBuuGDSk4cgSAPoF9uJh1kYuZF41CoNLA4UW0q+1Ch2AXFu6+SF5xWU0O5b4m5kQYv875EPeg2gwYP4Udy6PIyyym/ytNYPcGsteuxfWVlynu0JR3DrxDiFMI41qMg10fQG4iDJjPhdQi1p6I55n2gXg5WNf0kP4SRQgUFKoYnYMDUhoozMvFwdUaC0sV6QnGpSJd69YURkRgyM+nd2BvVEJlnBXYeUDDoRC+CopymNCnLun5JSzeG13Do7k/iTsVzoZPZ+DqF8CQCe+xZ1UMSTE59HquPg5ZF0me+SG2XbviMmoUUw9MJa8kj9mdZqNNioAji6H1CPBtwcdbz2OrteCVrrVrekh/iyIECgpVjI2j0WO0ICsToRI4e9mQnmBcKtK1aQ1lZRQcO4artSutPFvxe/TvxjhDbV6Ckjw4/i1N/Bzp39iLJXujSckpqsnh3HdcOXOK9R9Px8nLh2GT3+fghgRiT6fT5bFQ/L0l8WPHYenjg/dHs/kh8kf2JexjfMvxhNgHwa9jwc4Tuk8hLDaDHedSeKlLbRx1ljU9rL9FEQIFhSpG52B0FsrPygLAxdeWjKumGUGLFgitlrz9+wEYWHsg8XnxhKeGg09zCOgIh74EfSlv9g6lVG9g3k4lec3t4sqZU6ydNQ0Hdw8efucDTu5K4/yfibTqH0j9tm7Ev/oqsqAA3wVfcEmfxJywOXT27cxjdR+Dw19C0mno+xFSa8fsLedxs9PybIfAmh7WP6IIgYJCFXNNCAqyMwFw8balMLeUgpwSVFZW6Fq2JP/AnwD09O+JtYU1v1761Xhyh1chJx4i1hLoasMTbfz56egVLqUq0Veqm8sRFSLwyNSZRIXlcnxLHA06edOyfyBJ775H0clTeM36EBnoy8S9E7GztOP99u8jsi7D7plQpy/UG8Cu8ykcjc1kbI8QdJYWNT20f0QRAgWFKuba0lB+tmlG4GMDQHq88WFu07EjJZcuUZqYiE6jo4d/D7bEbqFYXwzBvcCtLvz5GUjJmB4hWFmo+GjL+ZoZzH3C5YiTrJtdIQLxkUXs/yWKWs3c6PxYKJnffEv2unW4jhqFfe/ezDs+r2KrqJUzbHoDENDvY/QSPtpygUAXHY+28qvpof0rFCFQUKhitDob1BYW5GeZZgTXdg6ZlodsO3YAKF8eGlB7ALkluey5sgdUKmg/BpIj4NIuXG21jOxSm61nkjkWl1EDo7n3MYrA+zh6ePLI1JmkJejZueIc3iGO9HquPgX795Hy8cfY9e6N66hX2Bu/t2KrqE9HOLseorZB98ng6Me6EwlcSM5lfO9QNOq74xF7d/RSQeEuQgiBztGJAtOMwNrOEp29ZfnOIcvgYCw8PMjffwCANp5tcLd2r1geavQw2HoaZwXAC52CcLPT8uGm80rymiom7nR4uQg8PGUGuZkqNi+KwMnLhn6vNEZ/OZaE18ejDQ3Fe9aHJBemMHn/ZEKdQo1bRQszYfNE8GoCrUdSWKLnk60XaOLrQP9GXjU9vH+NIgQKCtWAjYNjuRCA0WCcesUoBEIIbDp2IP/gQWRZGWqVmv61+rM/YT8ZRRlGB7O2L0H0H5B4Ep2lBeN6hhAWl8n2s8k1NKJ7j7jT4ay/JgJTZ1JSZMlvX5zE2lbDgDFNsCjO48orryC0WvwWfIHBypKJeydSrC/mky6foFVrYes7kJ8GAz4DtQVL9kWTlFPEOw/WR6Wq+RSU/5ZbEgIhhLMQYrsQIsr0fdNMC0IIfaWkNBsrlQcJIQ4LIS4KIX4yZTNTULjr0Tk4li8NAbj725FxNZ+yEmPYCNuOHTHk5FB4+jRgXB4qk2VsjtlsPKHFs2BpC39+DsCjLf2o5WbD7C3nKdMrAelulXIR8PLm4akzkdKKjZ+FIwQMfLUpOp2K+Ndeo+xqIr6ff47G25svT37J8ZTjTG03lUCHQLi4A8JXQsdx4N2U5JwivvzjEn0betIq0Lmmh/ifuNUZwVvATillCLDT9PtmFFZKSjOwUvlsYK6UMhjIBJ6/xf4oKNwR2Dg5k5dZsabvHmiPNEjSrhmM27UDIcqXh0KcQqjrXJeNl0zvSdaO0OIZiFgLWZexUKuY+EBdLqXm88MRJfTErRB3qpIITJmBSq3j189OUphXyoOjm+DooSN51mwKDh7Cc9o0dM2bcfDqQZacWsKQ4CE8WOtBKM6FX8eBax3oPAGAT7ddoMxg4K2+dWt4hP+dWxWCQcA3puNvMOYd/leY8hR3B67lMf5P5yso3MnYu7pTmJNNaUkxAO4BxmBjKXG5AKgdHbFq3Ii8ffvKzxlYeyBn089yIeOCsaDtyyBE+aygd30P2tZyZs72SLILSm/jaO4d4k6Fs/6jChHQaG35fcFJMpPz6TeyEe4B9mT++BOZq1bh/OyzOA4dQlphGpP2TaKWQy3eam16193xHmTHw6AFoLHi7NUcfjkWz9PtAglwsanRMf5/uFUh8JBSJpqOkwCPv2hnJYQIE0IcEkJce9i7AFlSymvBVOIBn7+6kRDiRdM1wlJTU2+x2woK1Yu9qxsAuWnGf6s2jlqs7S1Jicspb2PbpQtFp05RZvr3PKDWADQqDWui1hgbOPhCk+Fw/FvITUYIwdQHG5BdWMq8nZG3d0D3ADHhx8xEQGttx+avTpEck0Pv5xvgV9+Z/CNHSPrgA2w6d8L9jfHoDXre2vcW+aX5fNLlE3QaHcTuh6NfG4XarzVSSmZsOouDtYYx3e/cwHJ/xz8KgRBihxAi4iafQZXbSeN2hr/a0hAgpWwJPA7ME0L858AbUsrFUsqWUsqWbm5u//V0BYXbir2rOwA5qSmA0UDsEWBXPiMAsOvRA4Dc3bsBcLRypGdAT36L/o2iMlNYiY6vg76kfAdRfW97Hm3lz3cH47iYojiZ/VsuHj3E+o+m4+Tjy8NTZmBlY8e2pWe4ci6Tbv+rR+1m7pTExpIw5lUs/f3x+fRThFrN0oilHE48zKQ2kwh2CoaSAtg4BpwCofs7AOy+kMKBi+mM7RGCg05TswP9f/KPQmBKSt/wJp8NQLIQwgvA9J3yF9dIMH1HA38AzYB0wFEIcc3tzhdQcvQp3BPYu5mEIK1i9uoWYE9mUj4lRcZJsLZOHTQ+PuTt3FXe5qGQh8gtyWV73HZjgUtt43bSsGWQnw7A+N51sNaomfH72ds0mrubCwf38evcD/EIqs0jU2ZibWvP7u/OEx2eSsdHQqjX3ouyzEwujxwJKhV+X32J2s6OsKQwFoQvoF9QP4YEDzFebPcMY+rJgZ+DpQ2legMzfj9HLVcbnmwbULMDvQVudWloI/C06fhpYMP1DYQQTkIIrenYFegAnDXNIHYDD/3d+QoKdyO2zi4IlYrctIp3I/cAO5CQdsU4KxBCYNujO/kHD2Iw5S9o5dkKfzt/VkeurrhYp/FQWgiHFgDgaqtlTI9gdl9I5Y8LN333UjBxZs9Ofp//MV4hdRk2eTpaGxv2/RLF+UNJtB4QRJPufhhKSogfPYayxCR8FyzA0t+ftMI0Ju6diJ+dH1PbTUUIAVeOwqGFxh1dQZ0BWHEglkup+bzdr95d4zx2M26157OAXkKIKKCn6TdCiJZCiK9NbeoBYUKIkxgf/LOklNdeZSYCrwshLmK0GSy9xf4oKNwRqNRqbJ1dypeGANwDjLlqk6Ir7AR23XsgS0rIO2DcPSSEYFidYRxPOU50tikEtVso1B8EhxcbHZiAZ9oHEeiiY/pvZylVtpPelFM7trBl4Vz8GjZm2KRpaHU6jvwaw+nd8TTp6UfLfoFIKUl8ezKFx47hPetDdM2bUWYoY8LeCeSU5PBpl0+x0dhAST6sGwn2PtDrfQBScoqYvzOKbqFu9Kz/V+bRu4NbEgIpZbqUsoeUMsS0hJRhKg+TUr5gOv5TStlIStnE9L200vnRUsrWUspgKeXDUsriWxuOgsKdg72rO9mVhEBnb4mjh47ES9kVZS1boHJwIG/X7vKygbUHYiEsWBO5puJind+EklyjGACWFiom96/PpdR8Vh6Kq/7B3GUc37yR7Uu+IKhZS4ZMmIrGyooT2y8TtimWeh286DAsGCEEaV8sIOe333AbNw77fv0A+Oz4ZxxNOsrUdlMJdQ41XnD7u5BxCQZ/CVZGQZ+1+TwlZQamDmhQU8OsMu7euYyCwh2Ok5c3mYnmZi+vYAcSL2YhDcZ9FcLCAtsuncn74w9kmdF24GrtSjf/bmy8tLHCaOzZEEL7G5cmiowzip713OkU4sqc7ZGk5irvUNc4smE1u1csJrhVOwa9MRkLS0vO7r/Kn2suUru5O12fqIsQgqz160lbsACHoUNxGfkiADvidrD8zHIeDX2UAbUHGC94cSccXQJtR0FQJwDCYjNYeyKBEZ2DCHK9+7aLXo8iBAoK1YSTlw8F2VkUF+SXl3nVdqS4oIyMxIoyux490WdlURAWVl42PHQ4WcVZFZ7GAF0mQFGWUQwwLiO9N7ABRaV6Ptx0rvoHdIcjpeTPX1ax7/sV1O3QhQfHTURtoeH8oUR2rzqPfwNnej1nDP2Qf+QIiVOmomvTBq/33kUIQUx2DO8ceIdGro2Y0MroJEZhJmwYDa6h0GMKAHqDZOqGM3g5WDGqW3ANjrjqUIRAQaGacPb2BSDjanx5mXeIA4DZ8pBt504InY6cTRUP/VaerQh2DOb7899XBJrzbgr1BsKfX5TvIKrtZsvIzrVZeyKBQ9Hp1T2kOxYpJXtXLefg6h9o0LUnfUe/jtrCgsgjSez65hy+oU70HdkItYWK4ugY4se8iqWfH76fzUdYWlJQWsDrf7yOpcqSOV3nYKk2RbvZNAHyU2DoItAYcw5/f+QyZxNzmNy/3l2Ra+DfoAiBgkI14exjEoKECiGwd7VG52DJ1aiKgHQqa2vsuncnd+tWZKnRY1gIwRP1nuB8xnmOpxyvuGi3yVCaDwfmlheN6haMr5M1U9ZH3JeGY4NBz7ZFnxP261qa9O5Pn5GvolKpuXgshR0rzuEV7Ei/VxpjYammNCWKafFQAAAgAElEQVSFKyNGINRq/BZ9hdrBASkl0w5OIzo7mtmdZ+Np42m88Jl1cPpn6DIRvJsBkJ5XzCdbL9CulstdFV30n1CEQEGhmnBw90SlVpvZCYQQeNV2JPFilllb+3790Gdnk3/wYHlZ/1r9cdA6sOrcqoqG7nWh8XA4sgRyrgJgbalm2sAGRKXksWx/TPUO6g6jrLSU3+bNJmL3NtoOG06P515CqFREh6eyfekZPIPs6T+qMRpLNfrcXK68OJKyzEz8Fi3C0s+YNOb789+zKWYTo5uOpp13O+OFcxLht9fAp4XRqc/EjN/PUVBSxrRBDYxbSu8RFCFQUKgm1BYWOLh7ms0IAHzqOJKXWUxWSkF5mU3HDqjs7cn5fVN5mbWFNcNChrHr8i4S8xIrLtB1Ihj0sPfj8qIe9TzoVd+DeTuiuJpVWH2DuoMoKSpk3expRB3+k65PjaDDI08ihCD2dBpbl0TgFmDHg6ObYGllYfQVGPMqxRcv4jt/HtaNGgJw8OpBPj76Md38uvF8I1PMS4Me1r0IZcUwZBGojcs/+6PSWHsigZe61KaOh11NDbtaUIRAQaEacfbxIz3ePFqoXz1jiOIrZyuik6osLbHr1ZPcHTswFFfsABoeOhyJ5IcLP1RcwCnQGJn0+LdGL1cT7w6oj0Qy7dcz1TKWO4nC3BxWT3+HK2dO8cArr9GivzHizeUz6WxedBoXH1sGjGmCpbUF0mAg8a23KDh0CO8ZH2Dbybjz53LOZd7Y8wZBDkF82OlDVML0ODwwD2L2Qt+PwNUYO6ioVM/k9acJcrW5ZwzElVGEQEGhGnELCCIz8SqlRUXlZQ7u1ti7WnH5rHnqSft+/TDk55O3Z095mZetFz38e7A6cjV5JZViC3V+A1QaY8J0E75OOl7tEcLWM8lsPp3IvUpeRjo/vfcWKXHRDHz9bRp0McZsij+fwaavTuPkacPAsU3R6jRIKUmeNYucTZtxf2M8DoOMgpFXkseru15FCMFn3T8zOo2B0Xt41wxoMBSaPVl+z892RhGXXsCMwQ2x0qhv+5irG0UIFBSqEfegWkhpIPVyxdq9EAK/+i4kXMhEX1Zh3LVp0wYLNzey164zu8ZzDZ8jtySXXyJ/qSi084R2r8DpXyD+WHnxiE61aOBtz5QNZ8gqKKm+gdUQmUlX+WHqBHLSUhn61jSCW7UF4PLZdH5bcAoHN2sGjW2KlY0x+FvGsmVkfvsdTk/9D+fnjUs/Bmlg0r5JxObE8mmXT/GzMyWYL8yCNc+Bgw8MmGcMAQ6cT8ph8d5ohjX3pX2w6+0f9G1AEQIFhWrEI8gYaDclJtqs3L+eM6XFepJjKraRCgsLHAYPJm/vXkqTKzySG7o2pI1XG749+y3F+kqOYx1fA1sP2DoJTFtMNWoVHz3UmMyCEqb/dm/5FiTHXOLHqRMoKSrkkakz8W/YGIDY02lsWngaR3cdg19rhrWdcetn1rr1pHz8CXZ9H8DjrbfKjbtfnPiCP+L/YEKrCbTxamO8uJRG43B2AgxbBlbGbb56g2TS2tPYW2uY3L/e7R/0bUIRAgWFasTOxQ0rWztSYi+ZlfvUdUKoBJfPmC8POQ4bCgYD2RvM4y++0OgF0grTKhLcA2jtjKGQrxw2bnU00cDbgZe61GLN8fh7Jihd7Mnj/PTeW6g1Goa/NxvP2sa1++jwVDZ/dRpnbxszEcjZuo3EyZOxad8O79mzESrjo25LzBaWnF7CsJBhPFb3sYobnPgOzqyF7pPBr1V58dL90Zy4nMWUB+vhbHPvZtJVhEBBoRoRQuAeWIvkGHMh0Fpb4FXbgdjTaWblloGB6Fq2JHvNmgpHMqCNZxsauDRgecRy9AZ9xQlNnwDPRrDjXSitsEOM6R5CbTcbJq+LIK+4jLuZM3t2sm72NBw9PHl8+ie4+BqXci4dT2Hr4ghc/ewYNK4pVrbG5aC8fftJeOMNrJs0wfeLL1BZGh/g4SnhTN4/mebuzZncZnLF9s/EU7DpTQjqAh3Gld/3Ykoun2yLpFd9DwY3/cucWfcEihAoKFQzXiGhpMbFmBmMAWo1dSM9IZ+s5AKzcoeHhlESF0fhsYq1fyEELzR6gcu5l9kWt62isUoNvWdA1uXy0BMAVho1Hz3UhKvZhczafHcuEUkpObzuZ7YsnItvvYY8+t5sbJ1dAIgKS2br12dwD7QvNwwDFISFET9mDNrgYPwWfYVKpwOMO4TG7BqDp40n87rNQ6M2JZApzIKf/wfWzjBsqfHvCZTpDYz/+SQ2lmpmDml0T/kM3AxFCBQUqhmfug2QBgNXo86blQc1NRoeo8PNU6/a9+6NytaWzJ9/Nivv7t+dYMdgFoYvpMxQ6S2/VhcI7Qf7Pi13MgNoEeDEcx2CWHno8l23RGTQ69m5dCH7f/yWep26MXTSe2hND/ULh5PYvvQMXrUdGPBqE7TWxn3+hacjuDLyJTReXvh/vQS1vTFKaFZRFq/sfAWAhT0X4mTlZLqJAda9ZMw9/Mg3YFuR+XDR3mhOxmczfXBD3Oy0t3HkNYMiBAoK1Yx3nbogBAnnzTOK2btY4+Zvd4MQqHQ6HAYPJmfzlvJ8xgAqoWJU01HE5sTye/Tv5jfpMxMMZbBlklnxm31CqeNhy5urT5GRf3fsIiotLmLjnJmc3L6Z1oMeou+o11FbGN/gT/8Rz44VZ/Gu41TuLAZQfPEiV0aMQO3ggP/yZVi4GGcOxfpixu4eS2JeIp91/4wA+0pZxA7MhcjNxr+dX+vy4vNJOczbEUn/Rl482Nj79g28BlGEQEGhmtHqbHALCCLhwo2pJWs1cyM5Joe8TPMw0s5PPgFlZWT++JNZeQ//HtR3qc+XJ7+kVF9a6YQgo2/B2fUQtb282EqjZt6jzcguKGXS2lNmdoc7kfysTH55fzKXjh2h+3Mv0enxZxBCIKUkbFMMe3+MJLCRKw+OaoxGa1zGKY6OIe7ZZxEaDf4rlqPxNMYKMkgDU/ZP4XjKcWZ0nEEz92YVN4r+A3Z9AA0fgtYvlhcXleoZ92M4DtYapg9ueDuHXqPckhAIIZyFENuFEFGmb6ebtOkmhAiv9CkSQgw21a0QQsRUqmt6K/1RULhT8QmtT2LkefRl5obb2s2MyxEXjyWblVsGBmLbuTOZP/6IoaTiTV4IwZhmY0jIS2Bt1Frzm7R/FVxC4HdTaksT9b3tGd+7DlvPJPPLMfNwF3cSaZdj+f6d8aRejmXg+Ldp1udBAKRBcuCXixzeGENoW0/6jmyIhWUlEXj6KTBI/Jcvw9Lfv/x6847NY3PsZsY2H8sDQQ9U3CgzFlY/B651YMD8cn8BMCabOZ+Uy8cPNbmndwldz63OCN4CdkopQ4Cdpt9mSCl3SymbSimbAt2BAqCStYs3r9VLKcNvsT8KCnck/g0bU1pcRGKkuZ3AydMG90B7zh9MuuEcp6f+hz493Sz+EEAH7w40c2/G4lOLKSyrFFfIQgsPzoGsONj7idk5L3SqRdtazkzbeIa49HzuNKJPHOWHqW9iKCtj+LTZhLQyBn8z6A3s+vYcJ3ddoXF3X3o8VQ+VKTdwZREI+GYF2uCK0A/LIpaVJ5h5vuHzFTcqzoUfHjPGExr+PWhty6t2nE1mxZ+xPNchiG513W/PwO8QblUIBgHfmI6/AQb/Q/uHgM1SyoJ/aKegcE/h37ApKrWa6PCwG+rqtvUkPSGPVFNS+2vYtG+PNiSYjOXLkYYKD2QhBGObjyWlMIUVZ1aYXyyoMzR+FA7Mh+SKmENqleDTR5qiUgnG/HCC4jI9dwJSSo5v3sj62dNx9PTm8Zlz8KhlfKCXlerZsjiiPNF8x4dDECrj2/vficDaqLXMPTaXvoF9ebvN2xU7fgx6WPMCpF6Ah1eAS+3yc5Kyi3hz9Unqe9kzsW/obRv/ncKtCoGHlPJaUJMk4J8yOA8HfriubIYQ4pQQYq4Q4i/N80KIF4UQYUKIsNTU1L9qpqBwR6LV6fCp24DYEzcKQUgrD1QWgguHzGcFQghcRoygODKSvF27zOpaeLSgd0Bvlp1eRlL+dbOJPjONnrHrXoJKdgQfR2s+ebgJp+Kz+eAO8DrWl5Wxc+lCdq9YTO2WbRj+3mzsnI07qYryS/n1s5PEnEyj8/A6tOofVP5A/zsR2Bm3k2kHp9HBuwMzOs6oCCQHsPN9iNwCfWdD7W4V/TBIXvspnKJSA58/3gytxb0XS+if+EchEELsEEJE3OQzqHI7abRC/aUlSgjhBTQCtlYqngTUBVoBzsDEvzpfSrlYStlSStnSzc3tr5opKNyxBDVtQerlWHLTzZ3IrGw0BDVyJfJIEvpS88Qy9v36ofH3J3XhwhsMva+3fB2DNDD32Fyzcmxc4cG5kHQK9s0xq+rTwJMRnYL47lAcG09epaYozMtl7az3yncGDXx9EhorKwBy0gpZ+/ExkmKy6f18Axp19S0/r+hC5F+KwJHEI7y5900aujZkTtc5Fb4CAOE/GKOKtnwOWr1g1pc52y9wMDqdaYMaUNvNlvuRfxQCKWVPKWXDm3w2AMmmB/y1B/3fbVZ+BFgnpSx/RZFSJkojxcByoPVfnq2gcJdTq4Xxn3fUkT9vqGvQ2YfC3NIbjMbCwgLXkS9SfPacWVRSAB9bH55u8DSbYjYRnnKdea3+QGj0MOz9CBJPmlVNeKAuLQKcmLTmFJdS87jdJMdcYuVb40g4F0Gfl8cZdwaZQkCkxOWw+qNjFOSUMGhsU0JaVSwyFIaHE/fUUwiVmoBvvzETgePJxxm9azQB9gEs7LEQnUZXccNLu2HjGOOyWd+PzIzDW88ksWD3JR5r7ccjLf2qf/B3KLe6NLQReNp0/DSw4W/aPsZ1y0KVRERgtC9E3GJ/FBTuWFx8/HALCOL8gT031PnWdcLJU8ep3fE3vPk7DByIxtubtC8WmNkKwBiDyN3anQ8OfUCpodSsjr4fgc4F1r1sFn5Co1bxxePN0GrUvLzy2G0NQXFmz05+nPIm0mDg0Wmzadi1Z3ld7Kk01n16HAsLFUPfbIF3SMUmxPw//yTuuedROzgQsGoV2toV6/vhKeG8vONlPHQeLO61GAetQ8UNE0/CT08adwg9uhIqzRIupeYx/ueTNPF14L2BDap34Hc4tyoEs4BeQogooKfpN0KIlkKIr681EkIEAn7A9f8DVgkhTgOnAVfgg1vsj4LCHU3dDl1IjLpAdsqN9oBGXX1JicslOSbHvE6jwXXMGIoiIsjZvNmsTqfRManNJC5kXuDbM9+a30znDAO/gJQzsH2KWZWXgzWfDW/GpdR8xv14Ar2hev0L9GWl7Fj6JVsWzsU7tC5PzpqHV3CFUTZiTzybvjyFk6cNwya2wNnLprwud8cOrox8CUsfHwJWfoelb0Xcn5OpJ3lpx0u46dxY2mcpbrpKy8aZsbDqYbByhCdXl0cUBcgrLmPkd8fQWqj48skW96VdoDK3JARSynQpZQ8pZYhpCSnDVB4mpXyhUrtYKaWPlNJw3fndpZSNTEtNT0opb/88VUHhNhLazpgd6/yBvTfWtfXE0tqCE9sv31DnMHAA2nr1SP10jlkGM4CeAT3p4d+DL09+yeWc686t0xvajYYji+Gs+YS9Y4grUx+sz45zKXy01Xxba1WSm5HGT9MmcXLb77QaOIxhb09HZ28K86w3sOf7C+z5IRL/hi4Mfr0ZNg4Ve0ay1qwlfuw4tPXrEfDdt2jcK7Z1nk49zUvbX8LZypmlvZfirqu05TM/HVYOM6abfHIN2Fd4COsNkrE/nCA6NY/PH2uGt6N1tY39bkHxLFZQuI04uHvgW68hp3dvu2GZx9LKgsbdfIk+kUpavPk7kVCr8ZjwJqVXr5K5chXX83abt9GoNEw7OA2D+fsW9HjXmIR9wxjjW3IlnmoXwBNt/Fm0J5rV1eBsFn82gpVvjSPtchwDXnuLzk88i0ptfPsuzCvh1/nhROxNoFkvf/q93Lg8ZISUktTPvzCGkm7TmoBly1A7OpZf93jycUZuH4mj1pFlfZbhYVNpw2JhFqwcAllX4LEfwb2uWZ+m/3aWnedTmDao4T2baOa/ogiBgsJtpknvfmQnJxF76sSNdT380FipCdsUe0OdTbt22HTpTNqXX1KabG5Udte5M77leI4kHWHl2ZXmJ1pYwkPLjMe/PGPmdSyE4L2BDWhXy4W3157m4KX0Wx0eANJgIOy3dfw8/W20NrY8MWMOddp2LK9PT8jjlw/DSIrOoeez9Wk/LBiVyUdAlpSQOOlt0hYswGHIEPy++gqVTcVS0f6E/YzcPhIXaxeW9VmGp41nxY2LcmDlUEg5B8NXQUA7s34t2x/Dij9jGdEpiP+1DUDBiCIECgq3mZDW7dA5OBK+9bcb6qxsNDTu5sulEyk3OJgBeL79NrK0lOQPZtxQNyxkGN38ujHv+DzOZ1y31OMUCEO+hKsnYOOr5RnNwGg8/vLJ5vi76BjxbRgRCdncCjlpqayeMYU93y2ldos2PDFjTnkOAYBLJ1JY89Ex9GUGhoxvTmibige5PjeXyyNHkr1+Pa5jRuM1cwbCsiLUw5bYLYzZNYZAh0BWPLACL1uvihsX5xltAokn4eFvIKSXWb+2RCQx/fez9GngwaS+9262sf8PihAoKNxm1BYaGvfsS/Txo6Rejr2hvmlPf6x0Gg78EnXDDiLLgABcR48id/t2cnfsMKsTQjCt/TQctY5M2DvBPPwEQN3+xoxmp3+G/ea+B446S759rjX2VhY8s/wIsWn/PQyFlJIze3byzRujSLwYSa8XRzNw/Nvl4aP1ZQb2/xzFlkUROHnZ8MikVngE2ZefXxIXR+xjj1FwNAyvDz/EbdQoszwAqyNXM2HPBBq7NmZpn6W4WLtU3Lw4F75/FOKPGvMK1O1n1rd9Uam8+sMJmvg6Mu/RZuWzDwUjihAoKNQAzfsNxNLamkNrfryhzspGQ+sBQSREZt0QohrA5Zln0IaGkvT+dMoyM83qnKycmNlpJrHZsUw7OO3GaKOd3oCGw4xetufMZyTejtZ8+3wb9AbJ/5YdJinbPJHO31GQk83GT2eyZeFc3AKCeOqjz2nc44HyB3luRhHrPj3OyV1XaNTNl6FvNMfGscIonLdvHzEPP4I+NQ3/JYtxHFIRrcYgDcw/Pp9pB6fR3qc9X/X6CntL+0o3z4BvBsLlgzB0MTQwj3RzNDaDEd+GUcvNhm+ebY215f29Q+hmKEKgoFADWNva0eyBAUQePkDaTWYFDTp54+xtw4HVFykpMt/nLzQavD+cSVlmJomT37nhYd/Wqy2jmo7i9+jf+fbsdVtKhTBuKfVuBmueh9gDZtXB7raseLY1GXklDF98kKtZ180qbkLU0YOsGP8KMSeO0uXJ53jk3Zk4elQs98RFpPPTjCNkJObTZ0RDOj9aB7WF8dEjpSRtyRKuvDgSjZcXgWtWY9OuYl2/qKyIN/e8ydenv+ahOg/xeffPsbaotMsnNwmW9zPGVRq+Cho9ZNa3k1eyeG75UbwdrPnu+TY46DQo3IgiBAoKNUSL/oPRWuvY/e3XNzzMVWoVXR4LJTejiIPrLt1wrlX9+ni8MZ68XbvI/P77G+pfbPwivQJ6MefYHA4kmD/ssdTBE6vB0R9+GH6D53ETP0e+e6EN6XklPLr4IPGZN48RWVyQz5aFc9n4yQzsnF15ctZ8Wg4YiupauscSPft+iuS3L05i62jFI5NaEdyiYounPjeXhNdeJ/XTOdg90IfAH77H0rcinERaYRrPb3ue7XHbGd9iPFPbTkWjqvQgT7sIy/pA9hWjn0BoX7P+HYnJ4ImvD+Og07DyhTb3Raax/y+KECgo1BDWdva0f+RJLp8Ov2nYCe8QR5p08yNiTwJXzmfcUO/01FPYdOlMyqzZFBw334EkhOCDDh8Q7BjM+D3jOZt+XVIcGxf43zrQ2hv326deMKtu7u/EyhfakF1QyqOLDt1gM7gccZJv3hjN2X27aTtsOI/P+ARXv4pdOKmXc/l55lFO7Y6nUTdfHprYAkePirAPhadOETNkKLnbt+M2/nV85swpzy8MRm/hR399lMiMSOZ0ncMzDZ8xzxscsxe+7mE0ED+10Rg+ohJ7I1N5atlh3O21/PJSO8VX4B9QhEBBoQZp2rsfbv6B7F6xmKK8G/0p2w6uhaOHjh3Lz5Kfbe5IJoTAe9YsLLy8iB89mpL4BLN6nUbHgh4LsLe05+UdLxObHWt+cQdfeGo9CBUs7wuJp8yqm/g58v2IthSW6hmy8ADH4jIpLS5i14pF/DJ9MhZaLY9N/5gOjzxZnkpSrzdwbEssq2eFUVJYxoBXm9D50TrliWSkwUD60qXEPv4E0qAn4LvvcB0xovwhL6Vk1blVPLvlWbQWWlb2W0nPgJ7m/T7+HXw3BOw8YcRO8G1hVr0hPIEXvgmjlqstP49sh5eDIgL/hLjTU9fdjJYtW8qwsBvD+Soo3I0kXYrihylvENyqHQ+Om2j+5gukxeex5qMw3PzsGPRas/L19WsUR0cTO/wxNB7uBHz3nZnjFUBsdixPb3kaS7UlX/f+2jxvL0D6JaOxtSTXuGTkZx77MTYtn2eWH0GVGMWQwsOUZCTTvO9AOj72FBqtVXm7lLgcdn13nvT4PGo3d6frE6FY2VQs5ZTEx5P4zhQKDh3CrndvvKa/j9qhIuxDdnE27x98n21x2+jq15UZHWeYG4VLi2Dr2xC2FGp3N+YUqBQ2QkrJvB1RzN8ZResgZ5Y81RIHa8UmUBkhxDEpZcvry5UZgYJCDeNZO4T2jzxJ5KH9nNy26YZ6V19buj9Vj8RL2fyx8jzyurhA2lq18P1sPiWxcVx+7nn02eZ+AIEOgSzqtYjismKe3vw0FzLMl4FwqQ3PbQZrZ/hmAESsMat2lvmM4U/6J2wkJacQ+o6k81MjykWgpKiM/T9HsXpWGIW5JfQd2YgHXmxYLgLSYCBj1SqiBw6i6PRpPN+fhs/8eWYi8OfVPxm6YSi7Lu9ibPOxzO8231wEMqJhaS+jCLQfA4//YiYC+cVljPnhBPN3RvFwC19WPt9GEYH/gDIjUFC4AzAY9Kz/aDqxJ48zdOK7BDZtcUObo7/HcOTXGBp386XjIyE3zBzy9uwhfvQYtKGh+C1ZjIWTeQrx6KxoRmwfQWFZIXO7zqWNVxvzG+SlGiN1XjkEXSZS1HI0R39bx/HfN4AQNB8wjF9lHX4+kUyHYBfmPdKUjHNZHFp/ibysYhp29qHt4NporS3KL1l0IZKk6e9TGHYMm44d8Xp/Ghrvirg/eSV5zD8+nx8v/Egth1rM7DSTBi6VIoFKCeGrYMsk4xLWkK9uMAqfuZrNmO9PEJuez4QH6jKyc60b/jYKRv5qRqAIgYLCHUJJYQE/vjuRrOQkhk58F9/6Dc3qpZQcWH2Rkzuv0LCzD52G17nBMSp3924Sxo7DwtMTv6++QlsryKw+IS+BV3a8QlzO/7V399FRlXcCx78/MpNMSIZkQhIIJAECAuFNobxZ2ForIFqVlV00llWO267n9Gh112O7eGxX9+y2a9vT9dSu2lJaV3r2wPpS1F3tBm0BSxUM2PCS8BZDYpIGE/M+TOYlM7/9496YCQUTlGQC9/mcM2fufe6TzDO/c+/85j7PnfvU8tDCh1hfvL7/h2ZPiMj2+zn41u/Z11ZEMGLdMfUvvrKBMdnWFT/Pl9XxzIsVXBNwkxMRcid5WX7bdPKm9n1Dj7a30/yT/6Bt61aSvF5yv/VNMtau7TcWUFpbyg/f/SHN3c2sL17PAwsewOPq62qiox7+5wGoehMKPw9rf2Zd6dT7GjHlubdrePz/juEb7ebHJfNZWhT3IzPjz5hEYBiXAH9rCy/8yyN0ftTMzQ9upGj+on7bVZV3tr/PH3d8wOR52ay4e1a/b+BgTeBSd+99aCRC3nf/lTEr+99qwR/28/Ceh9lVt4uVk1by7aXfJsuTRSgQ4OAbr3PgtZcJdLQzOb2D5fkfMe62f4NZa0CE06c62P96DbWHW+h2we+SwxR9LodHb5lDjjeFWCBA29attPx8M9HOTnwlJeTc/41+4xaVLZU8ceAJ9jbupTirmO8s/Q5zc+b2NTAShH3PwFs/Ao3Cin+2ZhUb1deTfex0JxtfOkx5XTvXzczlB389j7Hp5vLQgZhEYBiXiEBHOy9+759orj3FsnXrWXzruo+vze91aGc9e144iTcrhZVfnc34KRn9tofrG2i4/36ClZVk/NVaxm3cSJLX+/H2mMZ49sizPFX+FNnRMXwleA3+fccIBc4wad58ltx6GwU5yfDS19DGQ/wp5y7eC93BByeDpKS5uGpFIbO/OJHNb9fw5O9O4tUIj8pJZu7cTqy1lbRly8j91jfxzOibc6C6o5qny5+mtKaUjJQMvn7l1ymZUUJS73uLRaFiu/Wr5/ZamHGjNf9yVt9ZTVNXkCd/e5Jt79YxJtXNozfP4pYrJ5iuoEEyicAwLiGRYJAdm37CsT/sJm/aDFb83b3kTi7qV6fx/Q52/OII/rYQs5dPYMmaIlLT+27QpuEwzU89TcumTSRlZpL9jfvwrVuHuN3EYlFqD/6Rt3/zIo2HjoAq7ZPcfGnd3Vy78BZEhHB3D8f3/okjpRW0tqeQOqqDq674gDklXyY5z5omMlxXR83m/yTw8sukhAJUjJ9O+M6vcdP61YzxuFFVyk6XsaVyC7vrd5PqSuWuWXexYfYGvMl2YgoH4PAL1pzCrdWQOxuu/26/CeYb2rvZ8nYNW96pJRKNccfiQv5h5XSy0pIxBm9IEoGIrAMeA4qBxap6zk9nEVkN/BhIAjarau9MZlOAbcBY4ABwp6qGB3pdkwgMJ1BVju7Zxa4tm+nu6mT6kmUsunkt46b2DRSHAhHe/d9THN7VgMs9ilnLJjD32olk5MT9eKuigqbHv49/fxntBRNonTOTuo5WugtbRisAAAfgSURBVP1djM7IZNY1X6Jpmouf1j5Hm7+Dz/esYkHXNYRrkolGlJxCL3OXepnmfxb3wWfpOROjMzSfzho33UdPgcvFmFWraLzuFp5oTOEPVR+Rnt7GjKlVtI/aR1OwgSxPFrfPuJ3bZ9xu3SwuFoX6/XDov+HwixDqgLwrrXshzbwJRo0iGImy+0Qzr5Q3UFrxIarKTfMm8ODK6UzOTjtf2IxPMFSJoBiIAT8DHjpXIhCRJOAEsBKoB8qAO1S1UkSeB36tqttE5KfAQVV9ZqDXNYnAcJJufxf7X32J8h2vEe7uxjchn2mLlpI/czY5k6eQnplF2+luDpTWUFXWRCymZBekMX6ykpLqJxRo5MPqozQer6Snp4ekaIzcrgCF2XmMX7iCUN4M2jWThtoAzbVdoMIZdwc12YdImxZkYUYm01tTyKppJbpvP6H3TwGQkhHBWxQj+doF1E2/kpMeD+XhFn7f+B7NwUZUhWhgKp7gQlZmL2JFbpBiqhnXeZjkmp0QaAGXB2atQef/DU1ZizjVEuBgXTtlNa3srW7FH+rBN9rNbQsLuPPqSeT7Rg8QLeOTDGnXkIjs4vyJ4GrgMVW93l5/2N70ONAMjFfVnrPrfRKTCAwnCp7xc2LvHo7t2U3D8Upi0SgASS4XqWMySHK70agS6PLTEwoAfce2uHJJGV1IStpk3O58QmeUSKxv3EFiUcb4PyArVM/Y4CnS/dXEOtpwxd3wLuyC4wWjqJ6axtHi0TRk9dAe7iJE34xo3miMxaEwS2JuvhiBzGAEV6gNt0Y+rtOsGbyjc9nrWsQB9+do7vHgD/UQ7un7P1Oy01halMUNc/K4eupY3EnmJ08Xw/kSgetclS+yiUBd3Ho9sASrO6hdVXviyidyHiJyD3APQGFh4fmqGcZly5OWzrzrVjPvutVEQkFOV52gpaGejqbTBP1dRCMRVBVPejqedC/pvhwUH7GYj1BACHRF0KiiKJ7RbjzpbtJHK6M760ltb0SbWom2hyGaBzqepMwMknxZhMb7qMmOUZMR5nSomUCog3yNUoDg8/jwpWQyyZXOtHCEfH8LSYEW69u+Yt3tNNWHevNodo3nxKipVPrTaDkTgVAPcyIxUpNHkZbsIt+XSuHYNIrzvOR6PQPGw7h4BkwEIvImMP4cmx5R1VfOUT4kVHUTsAmsM4Lhel3DGIncKR4KZs+jYPa8i/DfigasMRFY9hleQYBc+7F8gLrG8BswEajqioHqDKABKIhbz7fLWoBMEXHZZwW95YZhGMYwGo6OtzLgChGZIiLJQAnwqlqDEzuB3pkkNgDDdoZhGIZhWD5TIhCRW0WkHrgaeE1ESu3yCSLyOoD9bf8+oBQ4CjyvqhX2v/hH4EERqcIaM/jFZ2mPYRiGceHMD8oMwzAcwtyG2jAMwzgnkwgMwzAcziQCwzAMhzOJwDAMw+EuycFiEWkGaj/ln2cDH13E5lyuTJwGz8RqcEycBmco4zRJVXPOLrwkE8FnISL7zzVqbvRn4jR4JlaDY+I0OImIk+kaMgzDcDiTCAzDMBzOiYlgU6IbcIkwcRo8E6vBMXEanGGPk+PGCAzDMIz+nHhGYBiGYcQxicAwDMPhHJUIRGS1iBwXkSoR2Zjo9iSaiNSIyGERKReR/XZZloi8ISIn7WefXS4i8qQdu0MisiCxrR86IvJLEWkSkSNxZRccFxHZYNc/KSIbEvFehtp5YvWYiDTY+1W5iNwYt+1hO1bHReT6uPLL+tgUkQIR2SkilSJSISIP2OUjY79SVUc8gCTgfazpmJKBg8CsRLcrwTGpAbLPKvsBsNFe3gh8316+EfgN1mRTS4F9iW7/EMblC8AC4MinjQuQBVTbzz572Zfo9zZMsXoMaw7zs+vOso+7FGCKfTwmOeHYBPKABfayFzhhx2NE7FdOOiNYDFSparWqhoFtwJoEt2kkWgM8Zy8/B/xlXPkWtezFml0uLxENHGqq+hbQelbxhcbleuANVW1V1TbgDWD10Ld+eJ0nVuezBtimqiFVPQVUYR2Xl/2xqaqNqvqevdyFNTfLREbIfuWkRDARqItbr7fLnEyBHSJyQETuscvGqWqjvXwaGGcvOz1+FxoXp8frPrtL45e93R2YWAEgIpOB+cA+Rsh+5aREYPy55aq6ALgBuFdEvhC/Ua1zUXN98VlMXAb0DDAVuApoBH6U2OaMHCKSDrwE/L2qdsZvS+R+5aRE0AAUxK3n22WOpaoN9nMTsB3rFP3D3i4f+7nJru70+F1oXBwbL1X9UFWjqhoDfo61X4HDYyUibqwk8F+q+mu7eETsV05KBGXAFSIyRUSSgRLg1QS3KWFEJE1EvL3LwCrgCFZMeq9E2AC8Yi+/CtxlX82wFOiIO6V1gguNSymwSkR8dtfIKrvssnfW2NGtWPsVWLEqEZEUEZkCXAG8iwOOTRERrDnZj6rqv8dtGhn7VaJH04fzgTUSfwLrCoVHEt2eBMeiCOvqjINARW88gLHAb4GTwJtAll0uwFN27A4DCxP9HoYwNluxujQiWH2wX/00cQH+FmtAtAq4O9Hvaxhj9Ss7FofsD7S8uPqP2LE6DtwQV35ZH5vAcqxun0NAuf24caTsV+YWE4ZhGA7npK4hwzAM4xxMIjAMw3A4kwgMwzAcziQCwzAMhzOJwDAMw+FMIjAMw3A4kwgMwzAc7v8BlqHyAQWS4q0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as pt\n",
    "\n",
    "# Forward Euler on a hyperbolic problem: terrible idea. Oh well.\n",
    "\n",
    "dt = 0.2*h\n",
    "u = cl.clmath.sin(cl.array.to_device(queue, grid))\n",
    "for i in range(1800):\n",
    "    _, result_dict = precomp_knl(queue, u=u, h=h)\n",
    "    out = result_dict[\"out\"]\n",
    "    out[0] = out[-1] = 0\n",
    "    u = u + dt*out\n",
    "    \n",
    "    if i % 300 == 0:\n",
    "        pt.plot(u.get())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
